(Gelenbe Erol) Analysis and Synthesis of Computer
(Gelenbe Erol) Analysis and Synthesis of Computer
(Gelenbe Erol) Analysis and Synthesis of Computer
Synthesis of
Computer Systems
2ND EDITION
Published
Analysis and
Synthesis of
Computer Systems
2ND EDITION
Erol Gelenbe
Imperial College, UK
Isi Mitrani
University of Newcastle upon Tyne, UK
Distributed by
World Scientific Publishing Co. Pte. Ltd.
5 Toh Tuck Link, Singapore 596224
USA office: 27 Warren Street, Suite 401-402, Hackensack, NJ 07601
UK office: 57 Shelton Street, Covent Garden, London WC2H 9HE
For photocopying of material in this volume, please pay a copying fee through the Copyright
Clearance Center, Inc., 222 Rosewood Drive, Danvers, MA 01923, USA. In this case permission to
photocopy is not required from the publisher.
ISBN-13 978-1-84816-395-9
Printed in Singapore.
The book has been revised and extended, in order to reect important
developments in the eld of probabilistic modelling and performance
evaluation since the rst edition. Notable among these is the introduction
of queueing network models with positive and negative customers. A large
class of such models, together with their solutions and applications, is
described in Chapter 4. Another recent development concerns the solution
of models where the evolution of a queue is controlled by a Markovian
environment. These Markov-modulated queues occur in many dierent
contexts; their exact and approximate solution is the subject of Chapter 5.
Finally, the queue with a server of walking type described in Chapter 2 is
given a more general treatment in Chapter 10.
Erol Gelenbe
Isi Mitrani
February 2010
v
February 11, 2010 13:12 spi-b749 9in x 6in b749-fm
Contents
vii
February 11, 2010 13:12 spi-b749 9in x 6in b749-fm
Contents ix
Index 309
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
Chapter 1
Basic Tools of Probabilistic Modelling
1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
Fig. 1.1.
of a possible sample path for this process is shown in Fig. 1.1: customers
arrive at moments a1 , a2 , . . . and depart at moments d1 , d2 , . . .
An examination of the sample paths of a queueing process can disclose
some general relations between dierent quantities associated with a given
path. For instance, in the single-server system, if N (t1 ) = N (t2 ) for
some t1 < t2 , and there are k arrivals in the interval (t1 , t2 ), then
there are k departures in that interval. Since a sample path represents
a system in operation, relations of the above type are sometimes called
operational laws or operational identities (Buzen [1]). We shall derive
some operational identities in section 1.7. Because they apply to individual
sample paths, these identities are independent of any probabilistic assump-
tions governing the underlying stochastic process. Thus, the operational
approach to performance evaluation is free from the necessity to make such
assumptions. It is, however, tied to specic sample paths and hence to
specic runs of an existing system where measurements can be taken.
The probabilistic approach involves studying the stochastic process
which represents the system. The results of such a study necessarily depend
on the probabilistic assumptions governing the process. These results are
themselves probabilistic in nature and concern the population of all possible
sample paths. They are not associated with a particular run of an existing
system, or with any existing system at all. It is often desirable to evaluate
not only the expected performance of a system, but also the likely deviations
from that expected performance. Dealing with probability distributions
makes this possible, at least in principle.
We shall be concerned mainly with steady-state system behaviour
that is, with the characteristics of a process which has been running for
a long time and has settled down into a statistical equilibrium regime.
Long-run performance measures are important because they are stable;
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
P (S(t + y) = j|S(u); u t)
= P (S(t + y) = j|S(t)), t, y 0, j = 0, 1, . . . . (1.1)
The right-hand side of (1.1) may depend on t, y, j and the value of S(t). If,
in addition, it is independent of t, i.e. if
state i and of anything that happened before that time. This very important
property will be referred to as the memoryless property.
The probability pi,j (y), regarded as a function of y, is called the
transition probability function. The memoryless property immediately
implies the following set of functional equations:
pi,j (x + y) = pi,j (x)pk,j (y), x, y 0, i, j = 0, 1, . . . . (1.3)
k=0
These equations express simply the fact that, in order to move from state
i to state j in time x + y, the process has to be in some state k after time
x and then move to state j in time y (and the second transition does not
depend on i and x). They are the ChapmanKolmogorov equations of the
Markov process. Introducing the innite matrix P(y) of transition functions
pi,j (y), we can rewrite (1.3) as
That assumption, together with (1.3), ensures that pi,j (y) is continuous,
and has a continuous derivative, for all y 0; i, j = 0, 1, . . . (we state this
without proof).
A special role is played by the derivatives ai,j of the transition functions
at t = 0. By denition,
pi,i (y) 1
ai,i = lim , i = 0, 1, . . .
y0 y
(1.6)
pi,j (y)
ai,j = lim , i = j = 0, 1, . . . .
y0 y
Hence, if h is small,
In fact, since P(y) is a stochastic matrix (its rows sum up to 1), the rows
of P (y) must sum up to 0 for all y 0.
Let A = [ai,j ], i, j = 0, 1, . . . be the matrix of instantaneous transition
rates. Dierentiating (1.4) with respect to x and then letting x 0
yields a system of equations known as the ChapmanKolmogorov backward
dierential equations:
Either (1.10) or (1.11) can be solved for the transition probability functions,
subject to the initial conditions P(0) = I (the identity matrix) and P (0) =
A. In a purely formal way, treating P(y) as a numerically valued function
and A as a constant, (1.10) and (1.11) are satised by
H i (x)H
i (x + y) = H i (y), x, y 0. (1.14)
Any distribution function which satises (1.14) must fall into one of
the following three categories:
(i) H i (x) = 1 for all x 0. If this is the case, once the process enters
state i it remains there forever (properly speaking, the holding time
does not have a distribution function then). States of this type are
called absorbing.
(ii) H i (x) = 0 for all x 0. In this case the process bounces out of state
i as soon as it enters it. Such states are called instantaneous.
(iii) H i (x) is monotone decreasing from 1 to 0 on the interval [0, ) and
is dierentiable. States in this category are called stable.
From now on, we shall assume that all states are stable. Dierentiating
(x) = i H
Eq. (1.4) with respect to y and letting y 0 we obtain H i (x),
i
where i = Hi (0). Hence
i (x) = ei x ,
H x 0,
Hi (x) = 1 ei x , x 0. (1.15)
i = ai,i , i = 0, 1, . . . . (1.16)
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
From (1.15), (1.7) and the memoryless property it follows that the
probability that the process remains in state i for time x and then moves
to state j in the innitesimal interval (x, x + dx) is equal to
ei x ai,j dx, x 0, j = i.
Integrating this expression over all x 0 gives us the probability that the
next state to be entered will be state j:
ai,j ai,j
qi,j = ei x ai,j dx = = , i = j = 0, 1, . . . . (1.17)
0 i ai,i
We derived (1.15) and (1.17) under the assumption that the Markov
process was observed at some arbitrary, but xed, moment t. These results
continue to hold if, for example, the process is observed just after it enters
state i. Moreover, a stronger assertion can be made (we state it without
proof): given that the process has just entered state i, the time it spends
there and the state it enters next are mutually independent.
The behaviour of a Markov process can thus be described as follows:
at time t = 0 the process starts in some state, say i; it remains there for an
interval of time distributed exponentially with parameter i (average length
1/i ); the process then enters state j with probability qi,j , remains there for
an exponentially distributed interval with mean 1/j , enters state k with
probability qj,k , etc. The successive states visited by the process form a
Markov chain that is, the next state depends on the one immediately
before it, but not on all the previous ones and not on the number of moves
made so far. This Markov chain is said to be embedded in the Markov
process.
We shall conclude this section by examining a little more closely
the exponential distribution dened in (1.15). That distribution plays a
central role in most probabilistic models that are analytically tractable. It
owes its preeminent position to the memoryless property. If the duration
of a certain activity is distributed exponentially with parameter ,
and if that activity is observed at time x after its beginning, then the
remaining duration of the activity is independent of x and is also distributed
exponentially with parameter :
P ( > x + y) e(x+y)
P ( > x + y | > x) = = = ey = P ( > y).
P ( > x) ex
(1.18)
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
On the other hand, we have seen in the derivation of (1.15) that (excluding
the degenerate cases) the memoryless property implies the exponen-
tial distribution. There are, therefore, no other distributions with that
property.
Let 1 and 2 be two independent random variables with distribution
functions
F1 (x) = 1 e1 x ; F2 (x) = 1 e2 x ,
f1 (x) = 1 e1 x ; f2 (x) = 2 e2 x ,
f (x)dx = P ( = x) = P (min(1 , 2 ) = x)
= P (1 = x)P (2 x) + P (1 x)P (2 = x)
= f1 (x)dx[1 F2 (x)] + f2 (x)dx[1 F1 (x)]
= 1 e1 x e2 x dx + 2 e2 x e1 x dx
= (1 + 2 )e(1 +2 )x dx. (1.19)
The time until the rst completion is thus distributed exponentially with
parameter (1 + 2 ). The probability that activity 1 will complete rst is
given by
q1 = P (1 < 2 ) = f1 (x)[1 F2 (x)]dx = 1 /(1 + 2 ). (1.20)
0
Restriction (ii), together with (1.9) and (1.17), imply that the instan-
taneous transition matrix of the Poisson process has the form
0 0
0 0
A=
0
= (I + U). (1.22)
0
Here, I is the (innite) identity matrix and U is the matrix which has
ones on the rst upper diagonal and zeros everywhere else:
0 1 0 0
0 0 1 0
U=
0
.
0 0 1
(t)n n
P(t) = e(UI)t = et eUt = et U . (1.23)
n=0
n!
Now, the matrix Un has ones on the n-th upper diagonal and zeros
everywhere else. Therefore, the rst row of the matrix dened by the series
on the right-hand side of (1.23) is (1, t, (t)2 /2!, . . .). The probability of k
arrivals in the interval (0, t] is equal to
et (t)k
pk (t) = , k = 0, 1, . . . . (1.24)
k!
Because of the memoryless property, the probability of k arrivals in any
interval of length t is also given by (1.24). In a small interval of length h,
there is one arrival with probability p1 (h) = h+o(h). The probability that
there are two or more arrivals in an interval of length h is P>1 (h) = o(h).
These last properties (plus the memoryless one) are sometimes given as
dening axioms for the Poisson process.
Since the Poisson process is a Markov process, the holding times, i.e.
the intervals between consecutive arrivals, are independent and distributed
exponentially with parameter . This property too, can be taken as a
denition of the Poisson process; it implies the Markov property and
everything else. The expected length of the interarrival intervals is 1/.
Therefore, the average number of arrivals per unit time is . For that reason,
the parameter is called the rate of the Poisson process. The average
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
Fig. 1.2.
where we have used (1.24). We see that the processes resulting from the
decomposition are both Poisson (with rates 1 and 2 , respectively). Not
only that, these processes are independent of each other. This result, too,
generalises to arbitrary number of components.
The superposition and decomposition of Poisson processes are
illustrated in Fig. 1.2.
In analysing system performance, we frequently employ the technique
of tagging an incoming customer and following his progress through the
system. It is therefore important to know something about the system state
distribution that customers see when they arrive. In this respect, Poisson
arrivals have a very useful, and apparently unique property: they behave
like random observers. More precisely, let {S(t), t 0} be a stochastic
process representing the state of a queueing system. That system is fed with
customers by one or more arrival streams. Consider an arbitrary moment
t0 ; let S(t
0 ) be the system state just prior to t0 . Then, if the arrival streams
are Poisson, the random variable S(t 0 ) is independent of whether there is
an arrival at t0 or not (Strauch [8]). This is because S(t 0 ) is inuenced
only by the past history of the arrival processes, and that is independent of
whether there is an arrival at t0 (looking backwards in time, the interarrival
intervals are still exponentially distributed and hence memoryless).
Thus, an arrival from a Poisson stream sees the same system state
distribution as someone who just happens to look at the system, having
otherwise nothing to do with it (a random observer).
To appreciate this remarkable property better, let us take a contrasting
example where the arrival stream is decidedly not Poisson. Imagine a
conveyor belt bringing machine parts to an operator at intervals ranging
between 20 and 30 minutes; the operation performed on each part lasts
between 10 and 18 minutes. Two hours after starting the belt, a random
observer (the shop oor supervisor?) may well see the operator diligently
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
is the n-stage Erlang density function. The mean and variance of Tn are,
respectively, n/ and n/2 .
the limits
exist and are independent of the initial state, and (ii), these limits constitute
a probability distribution:
j = 1. (1.28)
j=0
P(y) = , y 0. (1.29)
In other words, if at any moment the process state has the steady-state
distribution, then it has the steady-state distribution at time y later, no
matter how large or small y is. The state distribution becomes invariant
with respect to time.
There are two important questions which arise in this connection.
First, under what conditions does a steady-state regime exist for a Markov
process? Second, how does one determine the steady-state distribution of
the process? We shall leave the question of existence until the end of this
section and concentrate now on the determination of the vector , assuming
that it exists.
Dierentiating (1.29) at y = 0, and remembering that P (0) = A, we
obtain a system of linear equations for :
A = 0. (1.30)
This is known as the system of balance equations, for reasons which will
become apparent shortly. Being homogeneous, that system determines the
vector up to a multiplicative constant; the normalising equation (1.28)
then completes the determination.
The balance equations have a strong intuitive appeal. To see this, let
us write the i-th equation in the form
ai,i i = aj,i j . (1.31)
j=0
j=i
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
Fig. 1.3.
0 0 = 1 1
(1 + 1 )1 = 0 0 + 2 2
(1.32)
(i + i )i = i1 i1 + i+1 i+1
(there are two arcs going out and two arcs coming into each cut, except
for node 0). Alternatively, cutting o the group of states (0, 1, . . . , i), for
i = 0, 1, . . . , we obtain an equivalent system of balance equations:
0 0 = 1 1
1 1 = 2 2
(1.33)
i i = i+1 i+1
(one arc going out and one arc coming into each cut). The general solution
of (1.33) is easily obtained by successive elimination:
0 1 . . . i1
i = 0 , i = 1, 2, . . . . (1.34)
1 2 . . . i
This leaves one unknown constant, 0 , which is determined from the
normalising condition (1.28):
1
0 0 1
0 = 1 + + + . (1.35)
1 1 2
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
Further, let Ri,j be the total average amount of time spent in state j, given
that S(0) = i:
Ri,j = lim Ri,j (t) = pi,j (u)du. (1.37)
t 0
(The rst inequality is obvious; the second follows from the Chapman
Kolmogorov equations (1.3); the irreducibility of the process ensures that
pr,i (v) > 0 and pj,k (w) > 0.) Hence, either all states are transient, or all
states are recurrent.
The case of all transient states can be disposed of quickly: if Ri,j is
nite for all i, j then, according to (1.37),
This important result is the point of departure for most analytic and
numerical studies of systems modelled by Markov processes.
Returning to the Birth and Death process considered earlier, we can
assert now that the necessary and sucient condition for existence of
steady-state is the convergence of the series appearing on the right-hand
side of (1.35); when it exists, the steady-state distribution is given by (1.34)
and (1.35). That assertion follows from the result above and from the fact
that the Birth and Death process is irreducible; the probability pi,j (t) of
moving from state i to state j in time t is obviously non-zero, for all i, j
and all t > 0.
Fig. 1.4.
i = i 0 , i = 0, 1, . . . . (1.41)
The necessary and sucient condition for the existence of a solution whose
elements sum up to 1, and hence for the existence of steady-state, is
< 1. When the system is in equilibrium, the number of customers in it is
distributed geometrically:
P (N = i) = i = i (1 ), i = 0, 1, . . . . (1.42)
The expectation, E[N ], and the variance, Var[N ], of that number are
given by
E[N ] = ii = /(1 ), (1.43)
i=1
and
Var[N ] = E[N 2 ] E 2 [N ] = i2 i E 2 [N ] = /(1 )2 . (1.44)
i=1
condition is, in both cases, < 2. We shall use the expected number
of customers in the system, E[N ], as a measure of performance. In the
M/M/1 system we have, from (1.43),
E[N ]M/M/1 = .
2
For the M/M/2 system, we nd rst
A similar inequality holds for any number of servers. The reason for the
worse performance of the M/M/c system is that its full service capacity
is not always utilised: when there are less than c customers in the system,
some servers are idle. The M/M/c system is, in its turn, more ecient than
c independent servers with separate queues (i.e. c M/M/1 systems), where
each new arrival joins any of the queues with equal probability. We leave
that comparison as an exercise to the reader. The lesson that emerges from
all this is that, other things being equal, a pooling of resources leads to
improved performance.
A limiting case of the M/M/c system is the system with innitely many
servers, M/M/. Clearly, there can be no queue of waiting customers here.
The solution of the balance equations is as in (1.48), top case, for all i:
i = (i /i!)0 , i = 0, 1, . . . . (1.50)
i = i 0 , i = 0, 1, . . . , K, (1.51)
Fig. 1.5.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
Our other example is the M/M/c/c system, where only customers who
nd idle servers are admitted. The classic application for such a model
is a telephone exchange with c lines. The steady-state distribution of the
number of busy servers is
i = (i /i!)0 , i = 0, 1, . . . , c, (1.52)
where
c
0 = (i /i!)1 .
i=0
T = E[tj+1 tj ].
Let Ti be the total expected amount of time that the process spends in
state i during a regeneration period (i = 0, 1, . . .). By the same argument
that led to equations (1.38) and (1.39) it can be shown that the long-run
fraction of time that the process spends in state i, and hence the steady-
state probability of state i, is given by
i = Ti /T, i = 0, 1, . . . . (1.54)
average time the process remains in state 0 on each visit is 1/, and the
average time it remains in state i (i > 0) on each visit is 1/( + ), we have
T0 = M0 /
(1.55)
Ti = Mi /( + ), i = 1, 2, . . . .
T0 = T1
(1.57)
( + )Ti = Ti1 + Ti+1 , i = 1, 2, . . . .
These last equations, together with T0 = 1/ (there is only one visit to state
0 during a regeneration interval) can be solved by successive elimination:
Ti = i /, i = 0, 1, . . . . (1.58)
Substituting (1.58) and (1.59) into (1.54) we nally obtain the desired
distribution
i = i (1 ), i = 0, 1, . . . .
Note that the above approach can be applied to the general Birth and
Death process as well, with obvious minor modications.
The rst rigorous proof of that relation was given by Little [6]; hence, it
is known as Littles result, or Littles theorem. However, the validity
of the result had been realised earlier and there were also proofs for some
special cases.
Consider an arbitrary queueing system in equilibrium, and let N, W
and be the average number of customers in the system, the average time
customers spend in the system and the average number of arrivals per unit
time, respectively. Littles theorem states that
N = W, (1.60)
Fix an arbitrary moment t in the steady-state. The customers who are in the
system at that moment are those who arrived before t and will depart after
t. Since the arrival process is stationary with rate and customers arrive
one at a time, the probability that there was an arrival at time t u is du.
Such an arrival is still in the system at time t with probability 1 Fw (u).
Therefore, point t u contributes an average of (1 Fw (u))du customers
to the ones present at time t. Integrating over all values of u yields
N= (1 Fw (u))du = W,
0
W = W0 + (N )(1/) + (1/).
Fig. 1.6.
W0 = m/2
= M2 /(2m). (1.66)
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
qj = j + j+1 + . . . ; j = 1, 2, . . . .
qj = qj1 / = qj1 ; j = 1, 2, . . . .
qj = j j = 1, 2, . . . , or j = qj qj+1 = j (1 ), j = 0, 1, . . . ;
T (n), the total amount of time the sample path remains at level n during
[a, b];
A(n), the number of jumps from n to n + 1 during [a, b] (i.e. the number
of arrivals who nd n customers in the system);
D(n), the number of jumps from n to n 1 during [a, b] (i.e. the number
of departures who leave n 1 customers behind). Clearly, A(n) > 0 for
n = m, m + 1, . . . , M 1 and D(n) > 0 for n = m + 1, m + 2, . . . M .
Note the similarity of form between (1.71) and the balance equations
(1.33) of the Birth and Death process. It should be realised, however,
that the content is very dierent. The relations (1.33) were between the
parameters i , i of a certain stochastic process, and the probabilities i ,
taken over the set of all sample paths of that process. Those relations could
be used to determine the probabilities. Here, on the other hand, we have
identities valid for any sample path of any queueing process. The equations
(1.71) can also be solved for p(n):
n1
(k)
p(n) = p(m) , n = m + 1, . . . , M, (1.72)
(k + 1)
k=m
M
p(n) = 1.
n=m
The fractions p(n) can thus be determined in terms of the fractions (n)
and (n). The latter are not, however, parameters of the process; they are
characteristics of the same sample path for which the former are sought.
Knowing the values of (n) and (n) for one sample path does not help to
nd the values of p(n) for another sample path.
Suppose now that the sample path N (t) is observed over longer and
longer periods of time, and that during those periods it attains wider and
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
exist and are non-zero for all n = 0, 1, . . . (except for 0 ). Continuing the
analogy with the Birth and Death process, one would naturally expect the
fractions pn to be the unique solution of the innite system of equations
pn = 1
n=0
(1.74)
n pn = n+1 pn+1 , n = 0, 1, . . . .
This is not necessarily the case, as can be seen from the following example.
Consider the sample path illustrated in Fig. 1.7. N (t) goes through
alternating busy and idle periods of unit length. During the i-th busy period
(i = 2, 3, . . .), it spends time /2n1 at level n, n = 1, 2, . . . , i 1, and the
rest of the time at level i (0 < < 12 ). It is easily seen that the limits (1.73)
for this sample path are
0 = 1, n = n = 2n1 /, n = 1, 2, . . .
n
p0 = 1/2, pn = /2 , n = 1, 2, . . . .
Fig. 1.7.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
If we were dealing with a Birth and Death process with the above
parameters, then a sample path should spend, in the long run, a fraction
1/(1 + 2) of its time in state 0, with probability 1. A pathological sample
path like the one in Fig. 1.7 may occur, but the probability of such an event
is zero.
T (1 + 1 + 21 + ) = T /(1 1 )
H = 1 + 2 + + j1 .
where we have used (1.75). The average number of class j customers in the
system is obtained, of course, from Littles theorem: Nj = j Wj .
It is intuitively clear that, with priority scheduling, higher-priority
customers receive better treatment at the expense of lower-priority ones.
The above expressions make that intuition quantitative. They also allow
one to address various optimisation problems. For instance, given the arrival
and service characteristics, and a cost function of the form
C= ci Wi ,
iR
completion, the customer with the shortest service time of those waiting
is selected and served to completion. Customers arrive in a Poisson stream
with rate ; the probability distribution function of their service times is
F (x).
This model can be reduced to the one with head-of-the-line priorities
by introducing an innity of articial customer classes, using the service
time x as a class index (for a rigorous derivation, service times should be
rst assumed discrete and then a limit taken). Customers of class x arrive
at rate x = dF (x); the rst and second moments of their service times
are, of course, x and x2 , respectively. The trac intensity for class x is
x = xdF (x). Substituting these parameters into (1.78) and replacing the
sums by integrals we obtain the conditional expected response time Wx of
a customer whose service time is x (Phipps [7]):
Wx = x
x
x+
2
+ u dF (u) 2 1 udF (u) 1 udF (u)
0 0 0
x
x+
= x + (M2 /2) 1 udF (u) 1 udF (u) ,
0 0
(1.79)
where M2 is the second moment of F (x) and x and x+ denote limits from
the left and from the right (if F (u) is continuous at point x, the two are
identical). The unconditional expected response time W is given by
W = Wx dF (x). (1.80)
0
Note the similarity between (1.81) and (1.78). The numerator in the
second term of (1.81) also represents expected residual service, this time
averaged over classes 1, 2, . . . , j only.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch01
References
1. Buzen, J. P. (1976). Fundamental operational laws of computer system
performance. Acta Informatica, 7, 167182.
2. Buzen, J. P. (1977). Operational Analysis: An Alternative to Stochastic
Modelling. Research Report, Harvard University.
3. Cinlar, E. (1954). Introduction to Stochastic Processes. Prentice-Hall,
Englewood Clis, New Jersey.
4. Cobham, A. (1954). Priority assignment in waiting-line problems. Operations
Research, 9, 383387.
5. Foster, F. G. (1972). Stochastic Processes Proc. IFORS Conference, Dublin.
6. Little, J. D. C. (1961). A proof for the queueing formula L = W . Operations
Research, 9, 383387.
7. Phipps, T. E. (1961). Machine repair as a priority waiting-line problem.
Operations Research, 9, 732742.
8. Strauch, R. E. (1970). When a queue looks the same to an arriving customer
as to an observer. Man. Sci., 17, 140141.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
Chapter 2
The Queue with Server of Walking Type
and Its Applications to Computer System
Modelling
2.1. Introduction
Several important classes of computer subsystems can be modelled in
a unied manner using a single server queue, whose service becomes
unavailable in a manner which depends on the queue length after each
service epoch. Such models are particularly useful in the study of the
performance of certain secondary memory devices (paging disks or drums,
for instance) and in evaluating the behaviour of multiplexed data commu-
nication channels.
In this chapter we shall rst examine the properties of the basic
theoretical model, and then develop various applications. This will provide
us with a more economical presentation of the results. The performance
measures of each application will thus be obtained as special instances of
the more general results which will be derived rst.
Section 2.2 will be devoted to the presentation and analysis of the
queue with server of walking type which serves as the metamodel for
the computer system models. We rst derive the stationary queue length
distribution related to a special Markov chain embedded in the general
queue length process. Then, using general results from Markov renewal
theory, we obtain the stationary probability distribution for the model at
arbitrary instants. We prove that the latter is identical to the stationary
distribution at instants of departure (and hence at the instants of arrival of
customers); this generalises a similar result (due to Khintchine) which has
been proved for the M/G/1 queue. We also show in this section how the
M/G/1 queues analysis can be immediately obtained from the preceding
43
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
associated with the queue length at those instants. Under the present
assumptions it is easy to see that {Qn }n0 is a Markov chain. Therefore
the pk , if they exist, must satisfy
p0 = p0
0 + p1 0
(2.1)
k
pk = p0
k + pkj+1 j , k1
j=0
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
or
The quantities U (x), V (x) are readily obtained. Notice that due to the
Poisson arrivals of rate we have
(y)k y
k =
e dS(y)
0 k!
so that
U (x) = E e(x1)S (2.3)
But
so that
E[S])
1 + (E[S]
lim G(x) = p0 .
x1 1 E[S]
Therefore
(1 E[S])
p0 = E[S])
1 + (E[S]
yielding
1 E[S] xE[e(x1)S ] E[e(x1)S ]
G(x) = E[S])
1 + (E[S] x E[e(x1)S ]
.
Notice that p0 > 0 implies E[S] < 1; this is the stability condition for a
queue with server of walking type.
Consider the stationary queue length distribution gk , k = 0, 1, . . . ,
measured at instants just after a departure occurs, i.e. when the server
has just left the service time block of Fig. 2.1. The following relation is
obtained because a departure takes place given that Q > 0 in Fig. 2.1:
k
gk = pkj+1 j /(1 p0 ), k = 0, 1, . . .
j=0
where
W (x) = k xk = E[e(x1)s ].
k=0
Let us now consider the single-server queue with Poisson arrivals of rate ,
and independent service times S of probability distribution function S(x).
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
Notice that if the queue length is zero after a departure, the following
departure will correspond to the rst customer which will arrive. Let F (x)
be the generating function, for |x| 1,
F (x) = pk xk .
k=0
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
Notice that G(x) of (2.5) is identical to F (x) when S = S.
The expressions (2.5) and (2.8) give the generating functions for
the stationary queue length probability distribution at instants of time
embedded in the queue length process for the queue with server of walking
type, and for the M/G/1 queue, respectively. What we really would like to
have for both systems is the stationary distribution dened as
which we shall derive using some general results from Markov renewal
theory.
P [Qn+1 = j, Tn+1 t | Q0 = i0 , Q1 = i1 , . . . , Qn = in ]
= P [Qn+1 = j, Tn+1 t | Qn = in ].
For our purposes the Qn will take integer values and the Tn will be real-
valued; furthermore, both will be non-negative and Q0 corresponds to the
initial state at time zero. In our queueing models Qn will be the queue
n
length at the instant 1 Ti . Notice that
which is the probability that at some instant t between two instants of the
Markov renewal process the state is j, given that it was i just after the
n
most recent instant. By the term instant we refer to a time 1 Ti , n 1.
The result we seek will be obtained by applying the key renewal theorem
[5]. It states that
v(j)
lim P [Qt = k] = A(j, k, t)dt (2.11)
t m 0
j
is the average time the process spends in state k between two successive
instants, given that it was in state j at the most recent instant. Also, m is
the average time between instants. Thus the right-hand side of (2.11) is
merely the average time spent in state k between each successive pair of
instants, divided by m.
Let us now apply this result to the two queueing systems which we are
examining here.
Stationary queue length process for the M/G/1 system: For this
system, we see that the queue length process just after each departure (see
Fig. 2.2) is Markov renewal. This is because the time between two successive
departures is totally determined by the state just after a departure, as is the
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
queue length just after the following departure. The stationary probability
v(j) in (2.11) is therefore replaced by pj of (2.7). Also we see from Fig. 2.2
that m of (2.12) is given by
We will now show that L(x) = F (x); i.e. the stationary queue length
distribution (2.11) is identical to the stationary distribution just after
departure instants, for the M/G/1 queue.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
Notice that
t
t xy
e dt (1 S(y))e x dy = xet(x1) (1 S(t))dt.
0 0 0
Therefore
L(x) = p0
1 + (x 1) e t(x1)
(1 S(t))dt
0
+ F (x) et(x1) (1 S(t))dt.
0
But
1 1
et(x1) (1 S(t))dt = + E[e(x1)S ]
0 (x 1) (x 1)
1
= [V (x) 1].
(x 1)
Hence
F (x)
L(x) = p0 V (x) + [V (x) 1] (2.17)
(x 1)
d
lim L(x) = kpk .
x1 dx
k=0
E[S](1 + CS2 )
lim E[Qt ] = kpk = E[S] 1 + (2.20)
t 2(1 E[S])
k=0
In fact, Khintchines result is that the stationary queue length distribution and the
stationary distribution at instants of arrival are identical; but the latter is identical to
the stationary queue length distribution at departure instants.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
Call qk the left-hand side of this expression and dene the generating
function
M (x) = qk xk , |x| 1. (2.24)
k=0
Consider this expression separately for each case of (2.22). Take rst
k j = 0; this contributes the following term to M (x):
k
p0 t (tx) p0 t(x1)
(1 S(t))e
dt = (1 S(t))e
dt
m 0 k! m 0
k=0
p0
= (U (x) 1).
m(x 1)
The two cases of (2.22) covered by k j1, j > 0, contribute the expression
1 (tx)kj
pj xj et (1 s(t)) dt
m j=1 0 (k j)!
k=j
1 j t
(tx)kj+1
+ pj x e [(1 S(t)) (1 s(t))] dt
mx j=1 0 (k j + 1)!
k=j1
G(x) p0 V (x) W (x)
= W (x) 1 + .
m(x 1) x
Therefore
p0 G(x) V (x) W (x)
M (x) = U (x) 1 + 1 W (x) 1 + .
m(x 1) p0 x
dP [S < t]
= P [S t]/E[S]
= [1 S(t)]/E[
S].
dt
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
For a proof the reader may see Cox [7]. Intuitively, S is the amount of
time that a person arriving at a bus-stop will have to wait if buses pass by
at epochs 0, S1 , S1 + S2 , . . . and all the Si are independent and identically
distributed with common probability distribution function S(t). This is also
the residual lifetime of Chapter 1, section 1.6.
Now notice that
k (t)k = E[e
(x1)S
]1 U (x) 1
x et [1 S(t)]/E[
S] = .
k=0 0 k! (x 1)E[S]
(x 1)E[S]
Therefore let Q denote the number in queue in stationary state for the
queue with server of walking type and let Q be the corresponding quantity
for the M/G/1 queue (2.28) leads immediately to the following important
identity:
where a(z) is the random variable representing the number of arrivals (from
the Poisson arrival stream of rate ) during an interval distributed as the
random variable z. Since (2.28) is a relationship between probability distri-
butions, (2.29) is an identity in the sense of the probability distributions.
We can now compute directly E[Q ] using the PollaczekKhintchine formula
(2.20) for E[Q]:
(E[S])2 (1 + KS2 )
E[Q ] = + E[s] + 0 t(1 S(t))dt.
(2.30)
2(1 E[S]) E[S]
Another deeper and more general result is concealed in (2.29). We shall
state this fact without proof; the result is due to GelenbeIasnogorodski
[11]. Let W be the limit as n of the waiting time Wn of the n-th
customer arriving at a queue with server of walking type and with general
independent interarrival times, and let W be the limit as n of Wn the
waiting time of the n-th customer arriving at the corresponding GI/G/1
queue. The service time of this GI/G/1 queue is S = s + T (as in the case
of the M/G/1 queue corresponding to the queue with Poisson arrivals and
server of walking type). The result obtained in [11] is that
W = W + S (2.31)
and that W and S are independent, the equality being in the sense of the
probability distributions of the random variables on the left- and right-hand
sides of (2.31).
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
Let us see how (2.31) implies the result given in (2.29) when arrivals are
Poisson. In this case, the number of arrivals to the system during disjoint
intervals of time are independent. As in (2.29), let a(z) denote the number
of arrivals in time z for the Poisson arrival process. Then (2.31) implies that
But a(W ) + a(s) and a(W ) + a(S) are the numbers of customers remaining
in queue just after the departure of a customer in stationary state for the
queue with server of walking type and for the M/G/1 queue, respectively.
We have shown that, for both systems, the stationary queue length
distribution is identical to the stationary queue length distribution just
after departure instants; therefore
Q = a(W ) + a(s)
Q = a(W ) + a(S),
and
Fig. 2.3. (a) The physical model. (b) The mathematical model.
type of Fig. 2.1. The instant at which the test is Q > 0? is performed
corresponds to the time when the beginning of the k-th sector passes under
the read/write heads. The service time s is Y /N , the time necessary for
transferring a page. The rest period T after a service is Y (N 1)/N , the
time necessary for the beginning of the k-th sector to return under the
read/write heads. Finally, the idle period S if the queue is empty after
the rest in Fig. 2.1 is simply the time Y for one complete rotation of the
PD. Therefore, the model of a PD sector queue will be a special case of the
queue with server of walking type with
where = Y (N + 1)/2N , since the total time for serving a page transfer
will be uniformly distributed over the set of values Y k/N, k = 1, 2, . . . , N ;
C 2 is the squared coecient of variation of this service time, so that
2 N
2 2
Y (N + 1) 2 Yk 1 2 Y
(1 + C ) = = (N + 1)(2N + 1)
2N N N 3 2N
k=1
and
2 2 2N + 1
1+C = .
3 N +1
/ = (N + 1)/2
so that the PD with sector queueing can support (N + 1)/2 times as much
page trac as the PD without sector queueing, provided the page trac is
uniformly distributed among all of the sectors.
These results are illustrated on Fig. 2.4, where we show the average
queue lengths n and n for a PD with eight sectors.
In fact, the sector queueing policy may be viewed as a shortest access
(or service) time rst scheduling policy; such policies tend to optimise the
performance of service systems, as we shall see in Chapter 6.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
Fig. 2.4. Average queue length for paging drum with and without sector queueing
(N = 8).
1 L, rm r() rM .
The clock time, or time necessary for moving one bit in the shift register,
may be varied at will by appropriate electronic circuitry between the two
limits, as indicated earlier. Therefore, if the need arises, a variable clock
rate can be implemented in this system.
To simplify the discussion, and with no loss of generality, we shall
assume that the k-th sector which we are examining begins at 1 = 1 and
ends at 2 = L/N . Notice that these are cell positions (each cell containing
one bit) rather than units of rotation time as was the case with the paging
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
drum. If we use the queue with server of walking type to model the k-th
queue we must take
L/N L L
s= r()d, T =
r()d, S = r()d
1 (L/N )+1 1
for the corresponding service time and idle times. The integrals in the above
expressions should be, in the strict sense, summations. The loss in accuracy
in treating as a continuous random variable will be insignicant, however,
since L/N can be expected to be of the order of magnitude of 100 bits.
In view of (2.31), and the corresponding results (2.29) and (2.30) in
the case of Poisson arrivals of transfer requests to each sector, we see that
system performance will be optimised by setting r() to its minimum value
rm leading to a minimisation of queue length and transfer times. If the
arrivals of transfer requests to the k-th sector form a Poisson stream of rate
k , we can use (2.30) to obtain its average queue length nk :
(k rm L)
nk = + k rm L/N + k rm L/2.
2(1 k rm L)
In the previous analysis we assumed that the charge-coupled device
speed could be varied only as a function of angular position but not of
queue length. The obvious conclusion was that its rotation speed should
be maintained at its highest possible level at all times. The performance
of this device can be improved, however, if its speed can be varied as a
function of angular position, and also of queue length. This possibility has
been analysed in [8].
Let us assume for the time being that transfers to and from the device
occur in blocks of L bits so that there is only one queue of transfer requests
(N = 1). Consider now an idle period for the device, that is one in which the
queue is empty. It is clear that during such periods the initial address = 1
should dwell as long as possible in the vicinity of the read/write head so
that as soon as a transfer request occurs it may advance at maximum speed
to the head in order to minimise the latency delay preceding the beginning
of the transfer. Furthermore, as soon as the initial address passes under
the read/write head, the information cells of the device should be moved
as quickly as possible to the vicinity of the read/write head once again if
no arrivals have occurred.
In order to examine this behaviour more closely, assume that an idle
period begins at an instant which we arbitrarily x at t = 0 just after the
cell = 1 passes under the read/write head. Let D(t) be the number of
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
cells separating the address at time t from the initial address, this number
being counted in the direction of the motion of the cells. Thus D(0+ ) = L.
We shall assume that, as long as there are no transfer requests, the address
= 1 visits the read/write head every T seconds (T being xed).
Suppose that f (t)dt is the probability that an arrival occurs, ending
the idle period, in the interval (t, t + dt) for some t 0. Then the average
distance, in number of cells to be traversed, to the starting address for the
arriving transfer request is
=
D D(t)f (t)dt.
0
Our problem is to choose the function D(t) which will minimise D, since as
soon as an arrival occurs the optimum strategy will be to rotate the charge-
coupled device using the minimum clock time rm . The average latency, or
delay before the arriving request can begin its transfer, is then rm D.
Since both f (t) and D(t) are non-negative quantities, D is minimised
simply by letting D(t) be as small as possible for each value of t. We know
that D(kT ) = 0, and D(kT + ) = L for k = 0, 1, . . . ; furthermore, D(t) is a
decreasing function for every other value of t. It cannot decrease any faster
than 1/rm because of the limitation on rotation speed, and any slower than
1/rM because of the need to refresh the contents of each cell at least each
rM time units. This leads immediately of the optimum form for D(t) shown
in Fig. 2.5. This form guarantees that D(t) is as small as possible within
the given constraints so that D is minimised.
The time at which the speed of rotation must be changed during each
rotation is easily obtained from
(T )/rM = L /rm
yielding
= rm (T LrM )/(rm rM )
( ) = (T LrM )/(rm rM ).
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
Fig. 2.5. Optimum form for D(t), the number of cells separating the initial address
from the read/write heads during an idle period of the charge-coupled device memory.
f (t) = et
when there are arrivals per second to the system. This yields
T
=
D ekT D(x)ex dx .
k=0 0
Therefore
= (L + 1) e 1 eT e
D +
(1 eT ) rm rM
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
T = LrM
which yields
= L/(1 eLrM ) 1/rM .
D
Fig. 2.6. A multiplied data communication channel with polling time y and xed
message (packet) transmission time Y .
The service time of Fig. 2.1 will take the value Y . As far as the k-th buer
queue is concerned, we can also identify two cases which yield the best and
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
These two limiting cases can be analysed exactly and will provide perfor-
mance bounds for the k-th buer queue.
2.4.1. Best and worst case analysis for the buer queues
Let us rst consider the best case analysis for the k-th buer queue. Here
we will have s = Y + y, T = (N 1)y and S = N y for the queue with
server of walking type model of the buer service mechanism. The average
queue length will be the measure of performance which we will examine;
let bk be this quantity for the best case. Using (2.30) we can write
[k (Y + N y)]2
bk = + k (Y + y) + k N y/2
2(1 k (Y + N y))
where k is the rate of arrival of packets to the k-th buer.
For the worst case we have again s = Y + y, while T = (N 1)(y + Y )
and S = N y + (N 1)Y . If we denote by Bk the worst case average queue
length we have, again using (2.30),
[k N (y + Y )]2
Bk = + k (Y + y) + k [N y + (N 1)Y ]/2.
2(1 k N (y + Y ))
We see that the dierence between the best and worst cases is due both
to the values of the arrival rate of packets which will saturate the system,
which are [Y + N y]1 and [N (y + Y )]1 , respectively, and to the additional
terms which appear in the formulae. On Fig. 2.7 we show the form of these
results.
Fig. 2.7. Range of values taken by the average queue length of the k-th buer.
E[S] = N y + (j + 1)Y
= N y + jY
E[S]
s = y+Y
2k [N y + (j + 1)Y ]2 1
bk (j) = + k [(N + 2)y + (j + 2)Y ].
2(1 k [N y + (j + 1)Y ]) 2
References
1. Adams, C., Gelenbe, D. and Vicard, J. (1977). An Experimentally Validated
Model of the Paging Drum. IRIA Research Report, No. 229.
2. Borovkov, A. A. (1976). Stochastic Processes in Queueing Theory.
Springer, New York.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch02
Chapter 3
Queueing Network Models
73
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
(i) the nature of each node: how many servers there are, how fast they are,
what scheduling strategy is employed there;
(ii) the nature of the jobs: their arrival patterns, their routing patterns, the
amounts of service they demand from nodes on their route. Typically,
there will be dierent classes of jobs in the system, with dierent
characteristics.
Fig. 3.1.
at this node represents the time users spend thinking between receiving a
response to one job and submitting the next.
Node 2 contains a single server representing the CPU. A queue may
form here and, this being a time-sharing system, let us say that the
scheduling discipline is processor-sharing (we shall discuss processor-sharing
strategies in detail later). Nodes 3 and 4 are also single-server nodes,
representing the drum and the disk, respectively. The scheduling strategy
at both nodes is FIFO (rst-in-rst-out).
To model the fact that job execution consists of alternating CPU and
Input/Output intervals (the latter correspond to either page or le record
transfers) we impose the following routing rules: after leaving nodes 1, 3 or
4 jobs always go to node 2; after leaving node 2 they may go to nodes 1,
3 or 4 with certain probabilities. These probabilities will be assumed xed
normally but may, in some applications, be allowed to depend on past job
history. Think times, CPU times and I/O intervals are also governed by
probabilistic assumptions.
There may be jobs of dierent types in the system. For example, k of
the M terminals may be reserved for users with short, I/O-bound jobs
while the others are occupied by users with long, CPU-bound jobs. This
can be modelled by introducing two jobs classes with dierent routing
probabilities, think and service time distributions.
What do we expect to learn from the model? Some system performance
measures one may be interested in are: the average response time (the
time between leaving node 1 and returning to it) for jobs of class i; the
proportions of time that the CPU, the drum and the disk spend servicing
jobs of class i, the marginal and joint distributions of queue sizes at the
various nodes, etc.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
We shall return to this model at the end of the chapter, after the
tools of analysis have been developed. We shall then be able to write
down expressions for these performance measures in terms of the system
parameters.
Case 1. There is a single job class. Jobs arrive into node i from outside in a
homogeneous Poisson stream with rate 0i ; i = 1, 2, . . . , N . The amount of
service they require at node i is distributed exponentially with mean 1/i .
After service at node i (i = 1, 2, . . . , N ) jobs take an exit arc to node j (j = 0
or j > i) with xed probability pij (pij 0, j pij = 1). Node i contains
ci identical servers of unit speed (the last is not an important restriction,
it is made only to avoid extra notation), with a common queue served in
FIFO order. All external arrival and service processes are independent.
The state of the QN at any given time is dened as the integer vector
n = (n1 , n2 , . . . , nN ), ni 0, i = 1, 2, . . . , N
Fig. 3.2.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
when it exists.
We begin by observing that node 1 is entirely unaected by nodes
2, 3, . . . , N . It behaves like the classic M/M/c queue with parameters =
01 , = 1 and c = c1 ; the stationary distribution of the latter exists when
< c and is given by (see Chapter 1)
p(n) = (n) (k); n = 0, 1, . . . (3.1)
k=0
where
k
k
(0) = 1, (k) = (j), (j) = min(j, c), j, k = 1, 2, . . . .
j=1
Thus (3.1) can be used to obtain the marginal stationary distribution p1 (n1 )
of the number of jobs at node 1.
If the stream of arrivals into node 2 is also Poisson, say with total rate
2 , then node 2 would also behave like an M/M/c queue with parameters
2 , 2 and c2 ; the stationary distribution of the number of jobs at node 2,
p2 (n2 ), would exist if 2 < c2 2 and we could again use (3.1) to write an
expression for it. Furthermore, if n1 and n2 were mutually independent in
the steady-state, we could obtain their joint distribution by multiplying
p1 (n1 ) and p2 (n2 ).
What constitutes the input into node 2? It is formed in general by
splitting o part of the output from node 1 (a fraction p12 ) and merging
it with the external arrivals into node 2. Since Poisson streams remain
Poisson after splitting and merging (if the streams are independent), it will
be sucient, and necessary, to show that the total departure stream from
node 1 must be Poisson in order for the total arrival stream into node 2 to
be Poisson.
Perhaps the best way to approach this problem is via the notion of
reversibility, introduced by Reich [20]. He observed that, in equilibrium, the
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
Going back to our feedforward network, the output theorem shows that
the total arrival stream into node 2 is Poisson with rate 2 = 02 + p12 01 .
Thus the marginal stationary distribution of the number of jobs at node 2,
p2 (n2 ) is given by (3.1) with = 2 , = 2 and c = c2 . The theorem
also shows that, at any time t in the steady-state, the number of jobs at
node 2 is independent of the number of jobs at node 1. This is because only
departures from node 1 prior to t inuence the state of node 2 at t and
those departures are independent of the state of node 1 at t.
These arguments carry through to all other nodes in the network.
The total arrival stream into node j is Poisson with rate
j1
j = 0j + i pij , j = 2, 3, . . . , N (3.2)
i=1
where i is the total arrival (and hence departure) rate for node i, i =
1, 2, . . . , j 1. The derivation of (3.2) is obvious; it takes into account the
external arrivals into node j and those parts of the departure streams from
other nodes which are directed to node j. Furthermore, at any moment in
the steady-state, the states of the various nodes are independent of each
other because we have shown that the past departure stream is independent
of the present state at a node. Therefore, the stationary distribution of the
network state is equal to
Fig. 3.3.
queue behind it will go to node 4 via node 3, in which case they will arrive
there before J and cause it to wait longer. Thus the conditional probability
of a long sojourn time at node 4 given a long sojourn time at node 1 is
higher than the corresponding unconditional probability, i.e. the two are
not independent.
There is only one known case of a QN where the sojourn times of a
given job at dierent nodes are all independent: N nodes strictly in tandem;
all except the rst and the last contain a single exponential server; the rst
node is an M/M/c queue and the last can be an M/G/c queue (Burke [5];
Reich [20]). Another curious aspect of this problem is that if waiting times
are dened to exclude service times, then even in the above case the waiting
times at dierent nodes are not independent (Burke [5]).
N
j = 0j + i pij , j = 1, 2, . . . , N. (3.4)
i=1
We saw a special case of the trac equations in (3.2); there, the associated
matrix was triangular and the equations always had a unique solution;
the steady-state distribution existed if, and only if, that solution satised
i > ci i (i = 1, 2, . . . , N ).
It is readily seen that if a general Jackson QN has a steady-state regime
then the corresponding trac equations have a solution. Indeed, they are
satised by the total rates of input (number of arrivals per unit time),
1 , 2 , . . . , N into nodes 1, 2, . . . , N . To justify that statement it is enough
to observe that, in the steady-state, i is also the total rate of output from
node i (i = 1, 2, . . . , N ). The right-hand side of (3.4) then contains the rate
of external input into node j (0j ), plus all the output rate fractions which
are directed to node j (i pij , i = 1, 2, . . . , N ) i.e. the total rate of input into
node j. Thus the existence of a solution to (3.4) is a necessary condition for
the existence of a steady-state distribution of the Jackson QN. A rigorous
proof of this can be found in [12].
Before examining the suciency of that condition we shall introduce a
classication of the individual nodes of the network. This follows loosely the
one adopted by Melamed [18]. A node is called open if any job which visits
it is certain (will do so with probability 1) to leave the network eventually.
A node is called closed if any job which visits it is certain to remain in
the network forever. A node is called recurrent if any job which visits it
is certain to return to it eventually (clearly all recurrent nodes are closed
but not vice versa). For example, in the network of Fig. 3.4, nodes 5 and
6 are open, nodes 3 and 4 are closed and recurrent, node 2 is closed and
non-recurrent (transient), and node 1 is neither open nor closed.
Let A be the set of open nodes in the network, B be the set of the non-
open nodes and R be the set of the recurrent nodes. It can be demonstrated
that the trac equations (3.4) have a solution if, and only if, 0j = 0 for
all j B [18].
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
Fig. 3.4.
We shall give an outline of the proof. Suppose that 0j > 0 for some
j B and that j is either recurrent or has a path leading from it to some
node r R (otherwise j would be open). In either case, a fraction (perhaps
all) of the external arrivals into j nd their way to some recurrent node
r and hence keep on visiting it ad infinitum. Therefore r saturates in the
long run; the trac through it does not balance and (3.4) does not have a
solution. If, on the other hand, 0j = 0 for all j B, one solution of (3.4)
can be obtained by setting j = 0 for all j B and solving only those
equations in (3.4) corresponding to j A. That will be possible because all
nodes j A are transient (in Markov chain terminology) and therefore the
submatrix of (3.4) associated with them has an inverse.
If the trac equations (3.4) have a solution, 1 , 2 , . . . , N , then,
necessarily j = 0 for all j B R [18]. This can be explained intuitively
by remarking that jobs may leave the set of nodes B R but never arrive
into it (from the denitions of A and R and from the fact that 0j = 0,
j B). Hence that set of nodes eventually drains of jobs completely and
the trac through it, when balanced, is zero.
Bearing in mind that pij = pji = 0 for all i A and j R, we can
summarise the above results in the following manner.
Theorem 3.1. The trac equations (3.4) have a solution if, and only if,
they are equivalent to the three independent sets of equations
j = 0j + i pij , j A, (3.5)
iA
j = 0, jBR (3.6)
j = i pij , j R. (3.7)
iR
Note that (3.5) always has a unique solution because its matrix has an
inverse (due to jA pij < 1 for all i A). (3.7), if present, has innitely
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
Corollary 3.1. The trac equations have a unique solution if, and only
if, all nodes of the network are open.
If the set R of the recurrent nodes is not empty, there is (after all the
jobs have drained from the nodes in B R) a constant number of jobs
circulating in it. Furthermore, R is split into non-intersecting equivalence
classes by the relation communicate (nodes i and j communicate if
there is a path from i to j and a path from j to i). There is a constant
number of jobs circulating in each of these communicating classes and the
system of equations (3.7) splits into independent subsystems, one for each
communicating class.
Thus, in order to understand the steady-state behaviour of general
Jackson networks it is important to study two special cases:
(i) open networks all of whose nodes are open (we shall call these networks
completely open);
(ii) closed networks consisting of a single communicating class, with a
xed number of jobs circulating inside (we shall call such networks
completely closed).
N
= p(n1 , . . . , nj 1, . . . , nN )I(nj >0) 0j
j=1
N
+ p(n1 , . . . , nj+1 , . . . , nN )j (nj + 1)pj0
j=1
N
N
+ p(n1 , . . . , ni + 1, . . . , nj 1, . . . , nN )i (ni + 1)
j=1 i=1
i=j
j < cj j , j = 1, 2, . . . , N (3.9)
then the steady-state distribution of the network state exists and has the
form
Proof. First we verify that (3.10) satises the balance equations (3.8). We
substitute (3.10) into (3.8) and use the identities (see (3.1)):
j (nj )
j (nj 1) = j (nj ), nj > 0;
j
j
j (nj + 1) = j (nj ), nj 0.
j (nj + 1)
The factors j (nj )/ k=0 j (k), j = 1, 2, . . . , N , cancel out and (3.8) is
reduced to
N
N
0j + j (nj )I(nj >0) (1 pjj )
j=1 j=1
N
N N N
j (nj ) j (nj )
= 0j I(nj >0) + j pj0 + I(nj >0) i pij .
j=1
j j=1 j=1
j i=1
i=j
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
This last equation always holds. Indeed, individual terms on the left- and
right-hand sides can be equated: from (3.4) we have
N
1
1 pjj = 0j + i pij , j = 1, 2, . . . , N,
j i=1
i=j
N
j (nj ) j (nj )
j (nj )I(nj >0) (1 pjj ) = 0j I(nj >0) + I(nj >0) i pij .
j j i=1
i=j
As an important aside, we should point out that the last two identities
mean, in eect, that
(i) the rate of transition out of state n, due to a job leaving node j, is
equal to the rate of transition into state n, due to a job arriving into
node j; and
(ii) the total rate of arrivals into the network is equal to the total rate of
departures from the network.
Property (i) is usually called local balance (to distinguish it from the
global balance equations (3.8)) and it appears to be intimately connected
with the existence of product-form solutions.
Having established that (3.10) satises (3.8), we next verify, by direct
summation, that it also satises the normalising equation
p(n) = 1
n0
when (3.9) holds. Jacksons theorem now follows from the theorem in
section 1.4 which states that if the balance equations of an irreducible
Markov process have a positive solution which satises the normalising
equation then the steady-state distribution of the Markov process exists and
is given by that solution. The state process of a completely open Jackson
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
QN is, indeed, irreducible; this follows from the fact that the state n = 0
is accessible from every state [18].
Jacksons theorem implies that the states nj of individual nodes
(j = 1, 2, . . . , N ) at a given moment in the steady-state are independent
random variables. This is an even more remarkable result than in the case
of feedforward networks because, a priori, it would seem that the nodes of
a general Jackson network have more opportunities to inuence each other.
Furthermore, as we shall see at the end of this section, the total input
process into a given node is no longer Poisson, in general. Yet the node
behaves as if it were!
Remark: The theorem states that conditions (3.9) are sucient for the
existence of a steady-state distribution. Clearly, they are also necessary
because if steady-state exists the total rate of output from node j is j (j =
1, 2, . . . , N ) and, since the servers there are occasionally idle, the rate of
output is less than cj j (which is what it would be if all the servers were
busy all the time).
Suppose now that we are dealing with a completely closed network
with K jobs circulating inside. The state process of the network is a nite
Markov chain, it is irreducible (since all nodes communicate) and therefore
always has a steady-state distribution. The form of that distribution was
discovered by Gordon and Newell [11] (although it can be derived as a
special case from one of Jacksons theorems).
The proof of this theorem is also by direct substitution of (3.11) into (3.8)
and verifying that the latter are satised.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
Since jobs are being served at node j for an average of 1/j and they arrive
there at rate j , the average number of jobs being served at node j is,
according to Littles theorem, j = j /j (j = 1, 2, . . . , N ). If there is only
one server at node j then j is its utilisation factor (the fraction of time the
server is busy). The total average number of jobs being served (not waiting
in queues) in the network is 1 + 2 + + N and hence, again according to
Littles theorem, the total average amount of service a job obtains during
its residence in the network is
N
E[S] = j .
j=1
ej = j /.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
So far we have not used the distribution of the network state; the same
arguments would apply, for example, if interarrival and service times had
general distributions. One needs the state distribution if one is interested
in the numbers of jobs at various nodes or the time jobs spend there. In
particular, the average number of jobs at node j is equal to
E[nj ] = kpj (k).
k=1
which gives, once more according to Littles theorem, the average response
time W (the time between the arrival of a job into, and its departure from
the network):
W = E[n]/.
Let us now look at some trac processes of jobs between nodes. Very
little is known about these and the results that are available are mostly of
a negative nature. For example, the total input process into a node is not,
in general, Poisson. To demonstrate this, consider the single-node network
of Fig. 3.5 (Burke [6]). There is a single server at node 1; upon completion
of service a job leaves with probability p10 and is fed back with probability
p11 = 1p10 . The trac equation is 1 = 01 +p11 1 , yielding 1 = 01 /p10 .
Steady-state exists when 1 < 1 and the system, as far as the queue size
distribution is concerned, is equivalent to an M/M/1 queue with trac
intensity = 1 /1 :
p(n) = n (1 ), n = 0, 1, . . . .
Now, the queue size distribution left behind by departing (not fed-back)
jobs is the same as that seen by jobs arriving from the outside; the latter
Fig. 3.5.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
The j-th term in the sum is equal to the probability that the j-th customer
in the queue will be the rst to be fed back, multiplied by the density
function of j service times (Erlang with parameters j, 1 ). Next,
g(x) = p(n)gn (x) = p11 1 e(1 01 )x
n=0
after substitution of p(n) and gn (x) and inverting the order of summation.
This gives
x
p11 1
G(x) = g(t)dt = [1 e(1 01 )x ].
0 1 01
We see that the mean of F (x) is 1/1 (as expected), but F (x) is not
exponential and hence the input stream is not Poisson.
This situation raises the question of what is the network state
distribution at the moments when jobs move from one node to another
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
(or arrive from the outside). If the total input into a node is Poisson, and
independent of the network state, then the network state distribution at
input instants is the same as the steady-state distribution. We saw in an
example that the input process is not necessarily Poisson but we also saw
that input jobs may still see the steady-state distribution. This is, in fact,
the case for any completely open Jackson network: jobs arriving into a node
(externally or internally) do not, in general, form a Poisson process but they
see the steady-state distribution of the network state (Sevcik and Mitrani
[23]). In a closed network with K jobs circulating in it, a job coming into a
node sees the steady-state distribution of a network with K 1 jobs. These
results are generalised in [23] to a large class of networks with many job
classes; the networks may be open with respect to some job classes and
closed with respect to others.
A generalisation of the output theorem holds for the processes of
departure from Jackson networks in equilibrium: the stream of jobs
leaving the network from node j is Poisson with rate j pj0 and its past
is independent of the network state. Moreover, these streams are mutually
independent [18].
(k1 , k2 , . . . , kR ) 0, (3.13)
where k = k1 + k2 + + kR .
It is not dicult to verify, by direct substitution, that the solution of
(3.13) which satises the normalising equation
p(k) = 1
k0
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
is given by
R
p(k1 , k2 , . . . , kR ) = (1 )k! (kr r /kr !), (k1 , k2 , . . . , kR ) 0 (3.14)
r=1
R
= p(r2 , . . . , rk )r1 I(r1 >0) + p(rj , r1 , . . . , rk )rj (3.15)
j=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
k
p(r1 , r2 , . . . , rk ) = (1 ) r i (3.16)
i=1
R
p(k1 , k2 , . . . , kR ) = (1 )k! (kr r /kr !) (3.17)
r=1
Fig. 3.6.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
Let f (x) be the probability density function of and let f (s) be its
Laplace transform. Since the Laplace transform of the node l service time
is l /(l + s), we can write
L
l
f (s) = Al bl [i /(i + s)]. (3.21)
l=1 i=1
where the increment h, the staircase steps di and their number m are chosen
so that F (x) approximates F (x) with the desired accuracy. Each of the
unit step functions I(tih) is the distribution function of a constant (ih)
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
where Arl = ar0 ar1 arl1 is the probability that a class r job reaches the
l-th stage of its service.
When the required service times are not distributed exponentially,
the stochastic process dened by the number (or vector of numbers) of
jobs in the system is not Markov and one cannot nd its steady-state
distribution by means of balance equations. However, if those distributions
are Coxian, the Markov property can be reinstated by a suitable rede-
nition of the system state. The new process can then be studied in the
usual way.
In the case of the processor-shared server, dene the system state
as a vector of vectors (v1 , v2 , . . . , vR ), where vr = (kr1 , kr2 , . . . , krLr )
is a vector whose l-th element is the number of class r jobs which are
in the l-th stage of their service. As dened, the system state forms a
Markov process because all stages are distributed exponentially. We can
therefore write a set of balance equations for the steady-state distribution
of (v1 , v2 , . . . , vR ). These equations take into account transitions out of
and into a state due to arrivals of class r jobs and due to completions
of stage l of a class r service (r = 1, 2, . . . , R; l = 1, 2, . . . , Lr ). The
solution of the balance equations, subject to the normalising equation, is
given by
R
Lr
kr krl
p(v1 , v2 , . . . , vR ) = (1 )k! r [(Arl /rl ) /krl !] (3.22)
r=1 l=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
= (r /r ) = r (Arl /rl ) .
r=1 r=1 l=1
R
R
= (1 )k! [(r /r )kr /kr !] = (1 )k! (kr r /kr !).
r=1 r=1
using the same notation as in the previous case (and assuming that < 1).
The product on the right-hand side is dened as 1 if (r1 , l1 ) = (0, 0).
Aggregating over all states such that the class index of the rst job in
the LCFS order is r1 , that of the second job is r2 , etc., gives
k
k
p(r1 , r2 , . . . , rk ) = (1 ) (rj /rj ) = (1 ) r j
j=1 j=1
Lr
R
p(v1 , v2 , . . . , vR ) = e [(r Arl /rl )krl (1/krl !)]. (3.24)
r=1 l=1
Steady-state exists for all values of the parameters. In this case, since
exp() factorises into a product of exp((r Arl /rl )) over all r and l,
the random variables krl (the number of class r jobs in the l-th stage of
their service) are mutually independent. Aggregation of (3.24) yields (3.19).
with probability
pir,0 = 1 pir,js (i, j = 1, 2, . . . , N ; r, s = 1, 2, . . . , R).
j,s
The pair (i, r) associated with a job at a node is called job state. The set of
job states is split into one or more non-intersecting subsets (or subchains)
in the following way: two job states belong to the same subchain if there is
a non-zero probability that a job will be in both job states during its life
in the network. Denote these subchains by E1 , E2 , . . . , Em (m 1). (For
example, if jobs never change class when they go from node to node, there
will be at least R subchains.)
It may be that some subchains are closed, having a constant number
of jobs in them at all times, while others are open with external arrivals
and departures. Moreover, the external arrival processes may be state-
dependent in a restricted way. Let S be the state of the network (to be
dened later), let M (S) be the total number of jobs in the network in state
S and let M (S, Ek ) be the number of jobs in subchain Ek when the network
is in state S. The external arrivals may be generated in either, but not both,
of the following two ways:
Type 1 node: The service requirements for all job classes are distributed
exponentially with mean 1/i . Jobs are served in order of arrival. The state
Si of the node is dened as the vector (r1 , r2 , . . . , rni ), where ni is the
number of jobs present and rj is the class index of the j-th job in the FCFS
order. There is a single server whose speed Ci (ni ) depends on the number
of jobs and satises Ci (1) = 1 (multiple servers can be modelled by setting
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
In this model, the equations play the role which the trac equations played
in Jackson networks. The quantity eir is proportional to the total arrival
rate of class r jobs into node i (i = 1, 2, . . . , N ; r = 1, 2, . . . , R). Since
pir,js = 0 when the job states (i, r) and (j, s) belong to dierent subchains,
there are in fact m independent subsystems.
ejs = p0,js + eir pir,js ; (j, s) Ek , k = 1, 2, . . . , m. (3.27)
(i,r)Ek
where:
otherwise, if there are external arrivals and they are of type (ii) then
m M(S,Ek )1
d(S) = k (n) ,
k=1 n=0
where G is the same constant as in (3.28); d(n) is dened in the same way
as d(S) in (3.28); the factor gi (ni ) depends on the type of node i (i =
1, 2, . . . , N ):
if node i is of type 1 then
R
n
i
nir
gi (ni ) = ni ! (eir /nir !) [i Ci (j)] ,
r=1 j=1
1/ir is the average required service time for class r jobs at node i (node
types 2, 3 and 4):
Lit
1/ir = (Airl /irl ).
l=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
where
(1 i )ni i if node i is of type 1, 2 or 4
pi (ni ) =
ei ni i /ni ! if node i is of type 3,
provided that i < 1 for nodes of type 1, 2 or 4. We see that in this case
the nodes behave like N independent M/M/1 (for types 1, 2 and 4) or
M/M/ (for type 3) queues.
Some remarks are in order concerning the assumptions, generality and
usefulness of the BCMP model. Clearly, the introduction of dierent job
classes and node types widens considerably the eld of application of the
model. We shall give two examples of systems which can be modelled as
BCMP but not as Jackson networks.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
N N
g(z) = gi (z) = ni i z ni (3.33)
i=1 i=1 ni =0
dened whenever the component series converge. g(z) will be called the
generating function of the network, and the factors gi (z) the generating
functions of the individual nodes (i = 1, 2, . . . , N ). Clearly, the coecient
of z K in g(z) is precisely our normalising constant G: that coecient, like G,
is the sum of terms of the type n1 1 n2 2 nNN, one term for each composition
n1 + n2 + + nN = K.
Denote by i (z) the partial products in (3.33):
which implies the following recurrence relation for the coecients Gi (j):
least one job at node 1. A similar statement is true, of course, for any other
node. Since we are dealing with geometric series,
gi (z) 1
g(z) = i zg(z)
gi (z)
Hence
j K
1
E[ni ] = GN (K j). (3.38)
GN (K) j=1 i
where i (0) = 1, i (j) = ji /[Ci (1)Ci (2) . . . Ci (j)], j 1 with the previous
notation for i . The network generating function is
N N
g(z) = gi (z) = i (j)z j
i=1 i=1 j=0
This recurrence relation, together with the initial conditions G1 (j) = 1 (j),
j = 0, 1, . . . , allow GN (K) to be computed in O(N K 2 ) steps.
To nd the utilisation of node i, Ui , we proceed as before: G Ui is the
coecient of z K in the series
gi (z) 1
hi (z) = g(z) = g(z)[1 di (z)] (3.39)
gi (z)
where di (z) is the inverse of gi (z). Denoting the coecients of hi (z) and
di (z) by Hi (j) and Di (j), respectively, the convolution (3.39) yields a
recurrence relation
j
Hi (j) = GN (j) GN (s)Di (j s). (3.40)
s=0
di (z)gi (z) = 1
which yields
j
Di (0) = 1, Di (s)i (j s) = 0, j 1.
s=0
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
or
j1
Di (j) = Di (s)i (j s), j 1. (3.41)
s=0
Thus (3.41) can be used to compute Di (j) and then (3.40) to compute
Hi (j). The utilisation of the i-th node is given by
and
K
E[ni ] = jpi (j), i = 1, 2, . . . , N.
j=1
i = ei GN (K 1)/GN (K)
The methods described so far generalise to networks with more than one
job class. The generating functions of such networks are multi-variate (there
is one variable for each job class if jobs do not change classes; one variable
for each subchain if they do). The normalisation constant and various
quantities of interest are obtained by multi-variate convolutions (Reiser [21];
Reiser and Kobayashi [22]; Wong [25]). Some results remain unchanged: for
example, the throughput ir of class r jobs through node i, in a closed
network with Kr jobs of class r circulating inside (r = 1, 2, . . . , R), is given
by
where Si is the state of node i, pi (Si ) is the probability of that state and
nir /ni is the fraction of server capacity allocated to class r jobs (for type 2
nodes), or the probability of a class r job being in service (type 1 or 4
nodes). Such a procedure would involve the computation of the normalising
constant and then of the marginal probabilities pi (Si ).
When we talk about response times in the context of a closed network,
we usually mean the time between leaving a certain node and returning
to it. For example, in a terminal-driven system the collection of terminals
is modelled by one node (of type server-per-job). Let that be node i and
suppose that there are Kr terminals of class r, r = 1, 2, . . . , R (in a heavily
loaded system, when the terminals are busy all the time, jobs can be
identied with terminals). The response time for a class r job is dened as
the interval between the job leaving its terminal (the user presses carriage
return) and returning to it (the keyboard unlocks). Denote the average
response time for class r jobs by Wir . Let ir be the throughput of class r
jobs at node i and let E[nir ] be the average number of class r jobs at node
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
i (users in think state). The average number of class r jobs in the rest of
the system is Kr E[nir ] and, by Littles theorem,
On the other hand, node i being of type 3, jobs do not wait there; the
average sojourn time for class r jobs is equal to their average service time
(or think time) 1/ir ; again by Littles theorem, E[nir ] = ir /ir . Hence
the particular job class) uniquely determine each other. Moreover, (3.48)
and (3.47), like (3.46), are valid under much more general assumptions than
those of the BCMP model.
Let us now take, as an example, the terminal system introduced at
the beginning of this chapter (see Fig. 3.1) and obtain for it expressions
for some performance measures of interest. The system consisted of M
terminals (modelled by a node of type 3), one CPU (a type 2 node), one
paging drum and one ling disk (type 1 nodes). Suppose that there is only
one job class and that on leaving the CPU jobs go to the terminals, the drum
and the disk with probabilities p1 , p3 and p4 , respectively, (p1 +p3 +p4 = 1).
On leaving the terminals, the drum and the disk, jobs go to the CPU with
probability 1. Let 1/i , i = 1, 2, 3, 4 be, respectively, the average think
times, the average CPU intervals, the average drum transfer times and
the average disk transfer times (the latter two include rotational and/or
seek delays). The corresponding distributions may be arbitrary Coxian for
i = 1, 2, but have to be assumed exponential for i = 3, 4 (see section 3.5).
The ow equations, (3.7) or (3.27), are
e1 = p 1 e 2 e 3 = p3 e 2
e2 = e1 + e3 + e4 e 4 = p4 e 2
U2 = 2 G4 (M 1)/G4 (M ),
2 = 2 U2 = G4 (M 1)/G4 (M ).
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
References
1. Barbour, A. D. (1976). Networks of queues and the method of stages. Adv.
Appl. Prob., 8(3), 584591.
2. Baskett, F., Chandy, K. M., Muntz, R. R. and Palacios, F. G. (1975). Open,
closed and mixed networks of queues with dierent classes of customers.
J.A.C.M., 22(2), 248260.
3. Baskett, F. and Palacios, F. G. (1972). Processor Sharing in a Central Server
Queueing Model of Multiprogramming with Applications. Proc. 6th Ann.
Princeton Conf. on Information Science and Systems, pp. 598603. Princeton,
New Jersey.
4. Burke, P. J. (1958). The output process of a stationary M/M/s queueing
system. Ann. of Math. Stat., 39, 1141152.
5. Burke, P. J. (1972). Output Processes and Tandem Queues. Proc. Symp.
Computer Communications Networks and Telecommunications, Brooklyn.
6. Burke, P. J. (1976). Proof of a conjecture on the interarrival-time distribution
in an M/M/1 queue with feedback. IEEE Trans. on Comm., 24(5), 175176.
7. Buzen, J. P. (1972). Queueing Network Models of Multiprogramming.
Ph.D. Thesis, Harvard University, Cambridge, Massachusetts.
8. Chandy, K. M. (1972). The Analysis and Solutions for General Queueing
Networks. Proc. 6th Ann. Princeton Conf. on Information Science and
Systems, pp. 224228. Princeton, New Jersey.
9. Chandy, K. M., Howard, J. H. and Towsley, D. F. (1977). Product form and
local balance in queueing networks. J.A.C.M., 24(2), 250263.
10. Cox, D. R. (1955). A use of complex probabilities in the theory of stochastic
processes. Proc., Cambridge Phil. Soc., 51, 313319.
11. Gordon, W. J. and Newell, G. F. (1967). Closed queueing systems with
exponential servers. Operations Research, 15, 254265.
12. Jackson, J. R. (1957). Networks of waiting lines. Operations Research,
15, 254265.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch03
Chapter 4
Queueing Networks with Multiple
Classes of Positive and Negative Customers
and Product Form Solution
4.1. Introduction
In papers dating from the end of the 1980s and early 1990s [3, 6], new
models of queueing networks were introduced, in which customers can be
either negative or positive. Positive customers are the ones that we
are used to when we model service systems: they enter a queue, wait and
then receive service, and then they move on to another queue and the
same thing may happen until they nally leave the network (or continually
cycling inside the network indenitely). However in this new model called
a Gelenbe Network or G-Network, a positive customer may mutate into
a negative customer when it enters another queue. A negative customer
vanishes if it arrives to an empty queue, and otherwise it reduces by one
the number of positive customers in the queue it enters. Furthermore,
negative customers do not receive service so that their only eect is to
reduce the amount of work at the queue which they enter or to destroy
other customers, hence the term negative.
It has been shown [6] that networks of queues with a single class
of positive and negative customers have a product form solution if the
external positive or negative customer arrivals are Poisson, the service
times of positive customers are exponential and independent, and if the
movement of customers between queues is Markovian. This chapter will
discuss the theory of G-networks as it applies to networks of queues with
multiple classes of positive and negative customers, with direct relations
of destruction among negative customers of certain classes and positive
customers of certain other classes. We will also allow changes among
customer classes, as is usual in such models. Of course, as indicated in
117
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
previous chapters of this book, the classical reference for multiple class
queueing network models is [2] and the related theory is discussed there, and
in other sources. Multiple class queueing networks which include negative
customers were rst developed in [19] and generalised in [20]. The extension
of the original model [6] to multiple classes has also been discussed in [12].
Some applications of G-networks are summarised in [18]. G-Networks
can be used to represent a variety of systems. The initial model [6] was
motivated by the analogy with neural networks [4,11]: each queue represents
a neuron, and customers represent excitation (positive) or inhibition (nega-
tive) signals. Indeed, signals in biophysical neurons, for instance in the brain
of mammals, also take the form of random trains of impulses of constant
size, just like customers travelling through a queueing network. Results sim-
ilar to the ones presented in this paper have been used in [9 and 25], where
signal classes correspond to dierent colours in images. Other applica-
tions, including to networking problems [17] have also been developed.
Another application is to multiple resource systems: positive customers
can be considered to be resource requests, while negative customers can
correspond to decisions to cancel such requests. G-Networks have been
applied to model systems where redundancy is used to protect the systems
operation against failures: work is scheduled on two dierent processors
and then cancelled at one of the two processors as soon as the work is
successfully completed at the other, as detailed in [8].
The single server queue with negative and positive customers has
been discussed in [7], while stability conditions for G-Networks were rst
obtained under general conditions in [10]. G-Networks with triggers which
are specic customers which can re-route other customers [14], and batch
removal of customers by negative customers, have been introduced in [15].
Additional primitives for these networks have also been introduced in [13].
The computation of numerical solutions to the non-linear trac equations,
which will be examined in detail below, has been discussed in [5].
In this chapter we focus on G-Networks with multiple classes of
positive customers and one or more classes of negative customers,
together three types of service centers and service disciplines:
Type 1: rst-in-rst-out (FIFO),
Type 2: processor sharing (PS),
Type 4: last-in-rst-out with preemptive resume priority (LIFO/PR).
With reference to the usual terminology related to the BCMP
theorem [2], we exclude from the present model the Type 3 service centers
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
The following constraint must hold for all stations i of Type 1 and classes
N R
of negative customers m such that j=1 l=1 P [j, i][l, m] > 0
for all classes of positive customers k and p, Ki,m,k = Ki,m,p . (3)
This constraint implies that a negative customer of some class m arriving
from the network does not distinguish between the positive customer
classes it will try to delete, and that it will treat them all in the same
manner.
For a Type 2 server, the probability that any one positive customer of
the queue is selected by the arriving negative customer is 1/c if c is the
total number of customers in the queue.
For Type 1 service centers, one may consider the following conditions
which are simpler than (2) and (3):
ik = ip
(4)
Ki,m,k = Ki,m,p
for all classes of positive customers k and p, and all classes of negative
customers m. Note however that these new conditions are more restrictive,
though they do imply that (2), (3) hold.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
i,k + +
i,k
qi,k = S (5)
i,k + m=1 Ki,m,k [i,m +
i,m ]
N
R
+
i,k = P + [j, i][l, k]j,l qj,l (6)
j=1 l=1
N
R
i,m = P [j, i][l, m]j,l qj,l (7)
j=1 l=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
where each gi (xi ) depends on the type of service center i. The gi (xi ) in (8)
have the following forms:
R
N R
N S
N
qi,k i,k (1 d[i, k]) = +
i,k +
i,m . (12)
i=1 k=1 i=1 k=1 i=1 m=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
Proof. Consider (6), then sum it for all the stations and all the classes and
exchange the order of summations in the right-hand side of the equation:
N R N R
N R
+
i,k = j,l qj,l P + [j, i][l, k] .
i=1 k=1 j=1 l=1 i=1 k=1
and,
R
N N
S
+
i,k +
i,m
i=1 k=1 i=1 m=1
R
N
N R S
N
+
= j,l qj,l P [j, i][l, k] + P [j, i][l, m] .
j=1 l=1 i=1 k=1 i=1 m=1
R
N S
N R
N
+
i,k +
i,m = j,l qj,l (1 d[j, l]).
i=1 k=1 i=1 m=1 j=1 l=1
The state dependent service rates for customers at service center j will
be denoted by Mj,l (xj ) where xj refers to the state of the service center
and l is the class of the customer concerned. From the denition of the
service rate j,l , we obtain for the three types of stations:
FIFO and LIFO/PR. Mj,l (xj ) = j,l 1{rj,1 =l} ,
x
PS. Mj,l (xj ) = j,l |xj,l
j|
.
Nj,l (xj ) is the deletion rate of class l positive customers due to external
arrivals of all the classes of negative customers
FIFO and LIFO/PR. Nj,l (xj ) = 1{rj,1 =l} Sm=1 Kj,m,l j,m
x S
PS. Nj,l (xj ) = |xj,l j| m=1 Kj,m,l j,m .
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
and
N
R
gj (xj + ej,l )
i,m = Mj,l (xj + ej,l ) P [j, i][l, m]. (17)
j=1 l=1
gj (xj )
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
Then for the three types of service centers, 1{|xi |>0} i (xi ) =
S
m=1 i,m 1{|xi |>0} .
Then, we substitute the values of Yi,m , Mi,k , Ni,k and Ai,k for a LIFO
station:
R
1{|xi |>0} i (xi ) = 1{|xi |>0} 1{ri,1 =k} (i,k + +
i,k )/qi,k
k=1
R
1{|xi |>0} 1{ri,1 =k} i,k
k=1
R
S
1{|xi |>0} 1{ri,1 =k} Ki,m,k i,m
k=1 m=1
S
R
+ 1{|xi |>0}
i,m 1{ri,1 =k} (1 Ki,m,k ).
m=1 k=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
We use the value of qi,k from equation (5) to obtain after some can-
cellations of terms:
R
1{|xi |>0} i (xi ) = 1{|xi |>0} 1{ri,1 =k}
k=1
S S
Ki,m,k
i,m +
i,m (1 Ki,m,k )
m=1 m=1
S
R
= 1{|xi |>0}
i,m 1{ri,1 =k}
m=1 k=1
R
and as 1{|xi |>0} k=1 1{ri,1 =k} = 1{|xi |>0} , we nally get the result:
S
1{|xi |>0} i (xi ) = 1{|xi |>0}
i,m . (18)
m=1
Similarly, we substitute the values of Yi,m , Mi,k , Ni,k , Ai,k and qi,k :
R
1{|xi |>0} i (xi ) = 1{|xi |>0} 1{ri, =k}
k=1
S S
i,k + Ki,m,k i,m + Ki,m,k
i,m
m=1 m=1
R
1{|xi|>0} 1{ri,1 =k} i,k 1{|xi |>0}
k=1
R
S
1{ri,1 =k} Ki,m,k i,m
k=1 m=1
S
R
+ 1{|xi|>0}
i,m 1{ri,1 =k} (1 Ki,m,k ).
m=1 k=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
We separate the last term into two parts, and regroup terms:
R
1{|xi |>0} i (xi ) = 1{|xi |>0} 1{ri, =k}
k=1
S S
i,k + Ki,m,k i,m + Ki,m,k
i,m
m=1 m=1
R
1{|xi |>0} 1{ri,1 =k}
k=1
S S
i,k + Ki,m,k i,m + Ki,m,k
i,m
m=1 m=1
S
R
+ 1{|xi |>0}
i,m 1{ri,1 =k} .
m=1 k=1
Conditions (2) and (3) imply that the following relation must hold:
R
S S
1{ri, =k} i,k + Ki,m,k i,m + Ki,m,k
i,m
k=1 m=1 m=1
R
S S
= 1{ri,1 =k} i,k + Ki,m,k i,m + Ki,m,k
i,m .
k=1 m=1 m=1
R
Thus, as 1{|xi |>0} k=1 1{ri,1 =k} = 1{|xi |>0} , we nally get the
expected result:
S
1{|xi |>0} i (xi ) = 1{|xi |>0}
i,m . (19)
m=1
R
gi (xi ei,k )
1{|xi |>0} i (xi ) = 1{|xi |>0} Ai,k (xi )(i,k + +
i,k )
gi (xi )
k=1
R
R
1{|xi |>0} Mi,k (xi ) 1{|xi |>0} Ni,k (xi )
k=1 k=1
S
+ 1{|xi |>0}
i,m Yi,m (xi ).
m=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
Finally we have:
R
S
xi,k
1{|xi |>0} i (xi ) = 1{|xi |>0}
i,m . (20)
|xi | m=1
k=1
R x
As 1{|xi |>0} k=1 |xi,k
i|
= 1{|xi |>0} , once again, we establish the relation we
need. This concludes the proof of Lemma 3.
Let us now turn to the proof of the Theorem 1. Consider the global
balance equation the networks considered is:
N R
(x) j,l + Mj,l (xj )1{|xj |>0} + Nj,l (xj )1{|xj |>0}
j=1 l=1
N
R
= (x ej,l )j,l Aj,l (xj )1{|xj |>0}
j=1 l=1
N
R
+ (x + ej,l )Nj,l (xj + ej,l )
j=1 l=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
N
R
+ (x + ej,l )Mj,l (xj + ej,l )d[j, l]
j=1 l=1
N
N R
R
+ Mj,l (xj + ej,l )(x ei,k + ej,l )
i=1 j=1 k=1 l=1
N
R
(j,l + Mj,l (xj )1{|xj |>0} + Nj,l (xj )1{|xj |>0} )
j=1 l=1
N R
gj (xj ej,l )
= j,l Aj,l (xj )1{|xj |>0}
j=1
gj (xj )
l=1
R
N S R
N
+ j,m Kj,m,l qj,l + j,l qj,l d[j, l]
j=1 l=1 m=1 j=1 l=1
N
N R
R
gi (xi ei,k )
+ j,l qj,l P + [j, i][l, k]Ai,k (xi ) 1{|xi |>0}
i=1 j=1 k=1 l=1
gi (xi )
N
N R
R
S
+ j,l qj,l P [j, i][l, m]Ki,m,k qi,k
i=1 j=1 k=1 l=1 m=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
N
N R
S
+ j,l qj,l P [j, i][l, m]Yi,m (xi )1{|xi |>0}
i=1 j=1 l=1 m=1
N
N R
S
+ j,l qj,l P [j, i][l, m]1{|xi |=0} .
i=1 j=1 l=1 m=1
After some substitution, we group the rst and the fourth terms of the right
side of the equation.
N
R
(j,l + Mj,l (xj )1{|xj |>0} + Nj,l (xj )1{|xj |>0} )
j=1 l=1
N
R
gj (xj ej,l )
= 1{|xj |>0} Aj,l (xj )(j,l + +
j,l )
j=1 l=1
gj (xj )
R
N S R
N
+ j,m Kj,m,l qj,l + j,l qj,l d[j, l]
j=1 l=1 m=1 j=1 l=1
R
N S S
N
+
i,m Ki,m,k qi,k +
i,m Yi,m (xi )1{|xi |>0}
i=1 k=1 m=1 i=1 m=1
N
S
+
i,m 1{|xi |=0} .
i=1 m=1
N R
We add to both sides the quantity j=1 l=1 j,l qj,l (1 d[j, l]) and
factorise three terms in the right side
N
R
(j,l + Mj,l (xj )1{|xj |>0} + Nj,l (xj )1{|xj |>0} ) + j,l qj,l (1 d[j, l])
j=1 l=1
N
R
gj (xj ej,l )
= 1{|xj |>0} Aj,l (xj )(j,l + +
j,l )
j=1 l=1
gj (xj )
R
N
S S
+ qj,l j,l + j,m Kj,m,l +
j,m Kj,m,l
j=1 l=1 m=1 m=1
S
N N
S
+
i,m Yi,m (xi )1{|xi |>0} +
i,m 1{|xi |=0} .
i=1 m=1 i=1 m=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
We substitute on the r.h.s, the value of qi,k in the second term. Then, we
cancel the term j,l which appears on both sides and we group terms to
obtain:
N
R
j,l qj,l (1 d[j, l])
j=1 l=1
R
N N
N
S
= +
j,l + 1{|xi |>0} i (xi ) +
i,m 1{|xi |=0} (21)
j=1 l=1 i=1 i=1 m=1
where
S
R
R
i (xi ) =
i,m Yi,m (xi ) Mi,k (xi ) Ni,k (xi )
m=1 k=1 k=1
R
gi (xi ei,k )
+ Ai,k (xi )(i,k + +
i,k ) .
gi (xi )
k=1
N
R
j,l qj,l (1 d[j, l])
j=1 l=1
R
N N
S
= +
j,l +
i,m (1{|xi |=0} + 1{|xi |>0} ).
j=1 l=1 i=1 m=1
As in the BCMP [2] theorem, we can also compute the steady state
distribution of the number of customers of each class in each queue. Let yi
be the vector whose elements are (yi,k ) the number of customers of class k
in station i. Let y be the vector of vectors (yi ). We omit the proof of the
following result.
Theorem 2. If the system of equations (5), (6) and (7) has a solution
then, the steady state distribution (y) is given by
N
(y) = hi (yi ) (22)
i=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
+ = + FP + + , = + FP + (24)
+ (I FP + ) = , (25)
+
= FP + . (26)
z = G(z) (29)
has a xed point z . This xed point will yield the solution of (25) and
(26) as:
(z ) = + z , + (z ) = (F (z )P + )n , (30)
n=0
and to notice that 0 Fi,k 1, and that (6), (7) now have taken the
form of the generalised trac equations (21). This completes the proof of
Proposition 2.
The above two propositions state that the trac equations always have
a solution. Of course, the product form (8) will only exist if the resulting
network is stable. The stability condition is summarised below and the proof
is identical to that of a similar result in [10].
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
4.5. Conclusion
In this chapter we have studied G-Networks. However, rather than develop
all of the theory, starting from networks with a single customer class,
we have dealt directly with G-Networks with multiple classes of positive
and negative customers. We have developed in detail both the existence
and uniqueness results for the steady-state solution of the model, and the
explicit product form solution. In the model considered, the service centers
are identical to the service centers considered in the BCMP theorem [2],
with the exception of the innite server case which is not considered.
However, all service times considered are exponentially distributed with
dierent service rates for dierent classes of positive customers.
Beyond this model, and the results discussed in [20] where multiple
classes of signals are allowed, where a signal is a generalisation of a negative
customer which has the ability to either destroy another customer or move
it to another queue, further extensions of these results can be expected to
emerge from future research.
We have mentioned applications of some of these results to algorithms
for colour texture generation [9, 25], using a neural network analogy where
colours are represented by customers of dierent types. The model we
have described, in a simpler single class version has also been applied to
texture recognition in medical images [21], and to optimisation problems in
computer-communication networks [22]. Other important characteristics of
these networks include their ability to approximate continuous and bounded
functions [23] which we think will lead to new developments in the eld of
stochastic networks and their applications.
References
1. Kemmeny, J. G. and Snell, J. L. (1965). Finite Markov Chains. Von Nostrand,
Princeton.
2. Baskett, F., Chandy, K., Muntz, R. R. and Palacios, F. G. (1975). Open,
closed and mixed networks of queues with dierent classes of customers.
Journal ACM, 22(2), 248260.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch04
Chapter 5
Markov-Modulated Queues
137
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch05
When the jumps of the random variable J are of size 1, i.e. when
jobs arrive and depart one at a time, the process is said to be of the
Quasi-Birth-and-Death type, or QBD (the term skip-free is also used
(Latouche et al., [7]). The state diagram for this common model, showing
some transitions out of state (i, j), is illustrated in Fig. 5.1.
The requirement that all transition rates cease to depend on the size of
the job queue beyond a certain threshold is not too restrictive. Note that
there is no limit on the magnitude of the threshold M , although it must be
pointed out that the larger M is, the greater the complexity of the solution.
Similarly, although jobs may arrive and/or depart in xed or variable (but
bounded) batches, the larger the batch size, the more complex the solution.
The object of the analysis of a Markov-modulated queue is to determine
the joint steady-state distribution of the environmental phase and the
number of jobs in the system:
That distribution exists for an irreducible Markov process if, and only if,
the corresponding set of balance equations has a positive solution that can
be normalised.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch05
Note that only transition (d) has a rate which depends on j, and that
dependency vanishes when j N .
Remark. The breakdown and repair processes could be generalised without
destroying the QBD nature of the process. For example, the servers could
break down and be repaired in batches, or a server breakdown could trigger
a job departure. The environmental state transitions can be arbitrary, as
long as the queue changes in steps of size 1.
In this example, as in all models where the environment state transi-
tions do not depend on the number of jobs present, the marginal distribution
of the number of operative servers can be determined without nding the
joint distribution rst. Moreover, since the servers break down and are
repaired independently of each other, that distribution is binomial:
i N i
N
pi, = ; i = 0, 1, . . . , N. (5.4)
i + +
Hence, the steady-state average number of operative servers is equal to
N
E(Xt ) = . (5.5)
+
The overall average arrival rate is equal to
N
= pi, i . (5.6)
i=0
This gives us an explicit condition for stability. The oered load must be
less than the processing capacity:
N
< . (5.7)
+
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch05
The only dependency on j comes from the fact that transitions (b), (c)
and (d) are not available when j = 0. In this example, the j-independency
threshold is M = 1. Note that the state (N, 0) is not reachable: node 1 may
be blocked only if there are jobs present.
N
1
F (x) = 1 i ei x .
i=0
back to the last established checkpoint; all transactions which arrived since
that moment either remain in the queue, if they have not been processed,
or return to it, in order to be processed again (it is assumed that repeated
service times are resampled independently).
This system can be modelled as an unbounded queue of (uncompleted)
transactions, which is modulated by an environment consisting of completed
transactions and checkpoints. More precisely, the two state variables, I(t)
and J(t), are the number of transactions that have completed service since
the last checkpoint, and the number of transactions present that have not
completed service (including those requiring re-processing), respectively.
The Markov-modulated queueing process X = {[I(t), J(t)]; t 0}, has
the following transitions out of state (i, j):
(a) to state (0, j + i), with rate ;
(b) to state (0, j)(i = N ), with rate ;
(c) to state (i, j + 1), with rate ;
(d) to state (i + 1, j 1)(0 i < N, j > 0), with rate ;
Because transitions (a), resulting from arrivals of faults, cause the queue
size to jump by more than 1, this is not a QBD process.
Aj = A; Bj = B; Cj = C, j M. (5.8)
Note that transitions (b) may represent a job arrival coinciding with
a change of phase. If arrivals are not accompanied by such changes, then
the matrices Bj and B are diagonal. Similarly, a transition of type (c) may
represent a job departure coinciding with a change of phase. Again, if such
coincidences do not occur, then the matrices Cj and C are diagonal.
By way of illustration, here are the transition rate matrices for the
model of the multiserver queue with breakdowns and repairs. In this case
the phase transitions are independent of the queue size, so the matrices Aj
are all equal:
0 N
0 (N 1)
..
Aj = A = 2 0 . .
.. ..
. .
N 0
Denoting
Also, let DjA , DjB and DjC be the diagonal matrices whose ith diagonal
element is equal to the ith row sum of Aj , Bj and Cj , respectively. Then
equations (5.9), for j = 0, 1, . . . , can be written as:
for j = M + 1, M + 2, . . . .
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch05
where Q0 = B, Q1 = A DA DB DC and Q2 = C.
Associated with equation (5.14) is the so-called characteristic matrix
polynomial, Q(x), dened as
Q(x) = Q0 + Q1 x + Q2 x2 . (5.15)
det[Q(xk )] = 0,
(5.16)
uk Q(xk ) = 0; k = 1, 2, . . . , d,
On the other hand, it is known that there cannot be more than d linearly
independent solutions (Gohberg et al., [4]). Therefore, any solution of (5.14)
can be expressed as a linear combination of the d solutions (5.17):
d
vj = k uk xjk ; j = M, M + 1, . . . , (5.18)
k=1
Proposition 5.1. The QBD process has a steady-state distribution if, and
only if, the number of eigenvalues of Q(x) strictly inside the unit disk, each
counted according to its multiplicity, is equal to the number of states of
the Markovian environment, N + 1. Then, assuming that the eigenvectors
of multiple eigenvalues are linearly independent, the spectral expansion
solution of (5.12) has the form
N
+1
vj = k uk xjk ; j = M, M + 1, . . . . (5.20)
k=1
1. Compute the eigenvalues of Q(x), xk , inside the unit disk, and the
corresponding left eigenvectors uk . If their number is other than N + 1,
stop; a steady-state distribution does not exist.
2. Solve the nite set of linear equations (5.11), for j = 0, 1, . . . , M , and
(5.13), with vM and vM+1 given by (5.20), to determine the constants
k and the vectors vj for j < M .
3. Use the obtained solution in order to determine various moments,
marginal probabilities, percentiles and other system performance mea-
sures that may be of interest.
u[Q0 + Q1 x + Q2 x2 ] = 0, (5.21)
by Q1
2 , it becomes
(a) Phase transitions leaving the queue unchanged: from state (i, j) to state
(k, j) (0 i, k N ; i = k), with rate aj (i, k);
(b) Transitions incrementing the queue by s: from state (i, j) to state
(k, j + s) (0 i, k N ; 1 s r1 ; r1 1), with rate bj,s (i, k);
(c) Transitions decrementing the queue by s: from state (i, j) to state
(k, j s) (0 i, k N ; 1 s r2 ; r2 1), with rate cj,s (i, k),
provided of course that the source and destination states are valid.
Obviously, if r1 = r2 = 1 then this is a Quasi-Birth-and-Death process.
Denote by Aj = [aj (i, k)], Bj,s = [bj,s (i, k)] and Cj,s = [cj,s (i, k)], the
transition rate matrices associated with (a), (b) and (c), respectively. There
is a threshold M , such that
Dening again the diagonal matrices DA , DBs and DCs , whose ith
diagonal element is equal to the ith row sum of A, Bs and Cs , respectively,
the balance equations for j > M + r1 can be written in a form analogous
to (5.12):
r1 r2
r1 r2
A Bs Cs
vj D + D + D = vjs Bs +vj A+ vj+s Cs . (5.28)
s=1 s=1 s=1 s=1
Similar equations, involving Aj , Bj,s and Cj,s , together with the corre-
sponding diagonal matrices, can be written for j M + r1 .
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch05
r1
r2
Q r1 = A D A DBs D Cs ,
s=1 s=1
where xk are the eigenvalues of Q(x) in the interior of the unit disk, uk are
the corresponding left eigenvectors, and k are constants (k = 1, 2, . . . , c).
These constants, together with the probability vectors vj for j < M , are
determined with the aid of the state-dependent balance equations and the
normalising equation.
There are now (M + r1 )(N + 1) so-far-unused balance equations (the
ones where j < M + r1 ), of which (M + r1 )(N + 1) 1 are linearly
independent, plus one normalising equation. The number of unknowns is
M (N + 1) + c (the vectors vj for j = 0, 1, . . . , M 1), plus the c constants
k . Hence, there is a unique solution when c = r1 (N + 1).
where H = Q Q1
r . Introducing the vectors y = x u, = 1, 2, . . . , r 1,
one obtains the equivalent linear form
0 H0
I 0 H1
[u, y1 , . . . , yr1 ] .. .. = x[u, y1 , . . . , yr1 ].
. .
I Hr1
As in the quadratic case, if Qr is singular then the linear form can be
achieved by an appropriate change of variable.
vj N +1 uN +1 j . (5.35)
This product form implies that when the queue is large, its size is
approximately independent of the environmental phase. The tail of the
marginal distribution of the queue size is approximately geometric:
vj = uN +1 j , (5.37)
eigenvalues that are near a given number. Here we are dealing with the
eigenvalue that is nearest to but strictly less than 1.
If (5.37) is applied to all vj , for j = 0, 1, . . . , then the approximation
depends on just one unknown constant, . Its value is determined by (5.13)
alone, and the expressions for vj become
uN +1
vj = (1 ) j ; j = 0, 1, . . . . (5.38)
(uN +1 1)
This last approximation avoids completely the need to solve a set of
linear equations. Hence, it also avoids all problems associated with ill-
conditioned matrices. Moreover, it scales well. The complexity of computing
and uN +1 grows roughly linearly with N when the matrices A, B and C
are sparse. The price paid for that convenience is that the balance equations
for j M are no longer satised.
Despite its apparent over-simplicity, the geometric approximation
(5.38) can be shown to be asymptotically exact when the oered load
increases.
lim vj = 0; j = 0, 1, . . . . (5.39)
1
where
The objective will be to show that both the exact distribution, where the
vectors vj are given by (5.20), and the geometric approximation, where
they are given by (5.38), have the same limiting distribution.
Consider rst the exact distribution. When all eigenvalues are simple,
the equations (5.20) and (5.39) imply that
lim k uk = 0; k = 1, 2, . . . N + 1. (5.43)
1
qG = 0; (q 1) = 1, (5.44)
Moreover, in view of (5.39), equation (5.45) holds if the lower index of the
summation is j = M (or any other non-negative integer), instead of j = 0.
Substituting (5.20) into (5.45) and changing the lower summation index
to j = M yields
N
+1
xM
k
lim k uk = q. (5.46)
1 1 xk
k=1
However, the rst N eigenvalues do not approach 1, while the last one,
xN +1 = , does. Hence, according to (5.43), the rst N terms in (5.46)
vanish and leave
N +1 uN +1
lim = q. (5.47)
1 1
Now, substituting (5.20) into (5.42), and arguing as for (5.47), we see
that only the term involving the dominant eigenvalue survives:
N
+1
h(s) = lim es(1)j k uk xjk
1
j=M k=1
N
+1
= lim k uk xjk es(1)j
1
k=1 j=M
N
+1
xMk e
s(1)M
= lim k uk
1
k=1
1 xk es(1)
N +1 uN +1
= lim . (5.48)
1 1 es(1)
1 1
h(s) = q lim =q . (5.49)
1 1 es(1) 1+s
The last limit follows from LHospitals rule. The Laplace transform
appearing in the right-hand side of (5.49) is that of the exponential
distribution with mean 1. Thus we have established the following rather
general result:
uN +1 1
= lim lim
1 (uN +1 1) 1 1 es(1)
1 uN +1
= lim , (5.50)
1 + s 1 (uN +1 1)
again using LHospitals rule.
The last limit in the right-hand side of (5.50) is simply the vector
q. This can be seen by arguing that the normalised left eigenvector of
the eigenvalue must approach the normalised left eigenvector of the
eigenvalue 1. Alternatively, multiply both sides of (5.47) by the column
vector 1:
N +1 (uN +1 1)
lim = 1. (5.51)
1 1
Hence rewrite (5.47) as
uN +1
lim = q. (5.52)
1 (uN +1 1)
Thus we have
1
h(s) =q = h(s). (5.53)
1+s
So, in heavy trac, the geometric approximation is asymptotically
exact, in the sense that it yields the same limiting normalised distribution
of environmental phase and queue size as the exact solution.
Fig. 5.4. Manufacturing blocking: Average node 1 queue size against arrival rate,
N = 10, = 1, = 1.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch05
at nodes 1 and 2 are the same. Hence, the busier node 1, the higher the
likelihood that the buer will ll up and cause blocking. Because of that,
the saturation point is not at = 1 (as it would be if node 1 was isolated),
but at approximately = 0.909.
The geometric approximation for the marginal distribution of the
environmental variable, I, indicating the number of jobs at node 2 and
whether or not node 1 is blocked, is given by (5.38) as q uN +1 /(uN +1 1).
Since there are two environmental states, I = N 1 and I = N , representing
N 1 jobs at node 2, the average length of the node 2 queue, L2 , is given by
N
1
L2 = iqi + (N 1)qN ,
i=1
where qi is the i+1st element of the vector q. Figure 5.5 compares the exact
value of L2 with that provided by the geometric approximation, for the same
parameters as in Fig. 5.4. It can be seen that this time the approximation
is relatively less accurate, and converges to the exact solution more slowly.
Intuitively, this is due to the fact that, in order to obtain an accurate value
for L2 , all elements of q need to be accurate. Whereas, in a heavily loaded
unbounded queue, only the tail of the distribution is important.
In Fig. 5.6, the average unbounded queue size is plotted against N .
Increasing the size of the nite buer enlarges the environmental state
Fig. 5.5. Manufacturing blocking: Average node 2 queue size against arrival rate,
N = 10, = 1, = 1.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch05
Fig. 5.6. Manufacturing blocking: Average node 1 queue size against N , = 0.8, = 1,
= 1.
Fig. 5.7. Breakdowns and repairs: Average queue size against arrival rate, N = 10,
= 1, = 0.05, = 0.1.
Fig. 5.8. Breakdowns and repairs: Average queue size against number of servers, = 6,
= 1, = 0.05, = 0.1.
5.11. Remarks
The presentation in this chapter is based on material from [8, 10, 11]. It
is perhaps worth mentioning that there are two other solution techniques
that can be used in the context of Markov-modulated queues. These are
the matrix-geometric method (Neuts, [12]) and the generating functions
method (as applied, for example, in [9]). However, we have chosen to
concentrate on the spectral expansion solution method because it is
versatile, readily implementable and ecient. A strong case can be made
for using it, whenever possible, in preference to the other methods [10].
An additional point in its favour is that it provides the basis for a simple
approximate solution.
The geometric approximation is valid for a large class of heavily loaded
systems. The arguments presented here do not rely on any particular model
structure. One could relax the QBD assumption and allow batch arrivals
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch05
References
1. Buzacott, J. A. and Shanthikumar, J. G. (1993). Stochastic Models of
Manufacturing Systems, Prentice-Hall.
2. Daigle, J. N. and Lucantoni, D. M. (1991). Queueing systems having phase-
dependent arrival and service rates, in Numerical Solutions of Markov
Chains, (ed. W. J. Stewart), Marcel Dekker.
3. Gail, H. R., Hantler, S. L. and Taylor, B. A. (1996). Spectral analysis of
M/G/1 and G/M/1 type Markov chains, Adv. in Appl. Prob., 28, 114165.
4. Gohberg, I., Lancaster, P. and Rodman, L. (1982). Matrix Polynomials,
Academic Press.
5. Jennings, A. (1977). Matrix Computations for Engineers and Scientists,
Wiley.
6. Konheim, A. G. and Reiser, M. (1976). A queueing model with finite waiting
room and blocking, JACM, 23(2), 328341.
7. Latouche, G., Jacobs, P. A. and Gaver, D. P. (1984). Finite Markov chain
models skip-free in one direction, Naval Res. Log. Quart., 31, 571588.
8. Mitrani, I. (2005). Approximate Solutions for Heavily Loaded Markov
Modulated Queues, Performance Evaluation, 62, 117131.
9. Mitrani, I. and Avi-Itzhak, B. (1968). A many-server queue with service
interruptions, Operations Research, 16(3), 628638.
10. Mitrani, I. and Chakka, R. (1995). Spectral expansion solution for a class
of Markov models: Application and comparison with the matrix-geometric
method, Performance Evaluation.
11. Mitrani, I. and Mitra, D. (1991). A spectral expansion method for random
walks on semi-infinite strips, IMACS Symposium on Iterative Methods in
Linear Algebra, Brussels.
12. Neuts, M. F. (1981). Matrix Geometric Solutions in Stochastic Models, John
Hopkins Press.
13. Neuts, M. F. and Lucantoni, D. M. (1979). A Markovian queue with N servers
subject to breakdowns and repairs, Management Science, 25, 849861.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
Chapter 6
Diusion Approximation Methods
for General Queueing Networks
6.1. Introduction
Although considerable progress has been made in obtaining exact solutions
for large classes of queueing network models, one particularly simple type of
network, an arbitrary network with rst-come-rst-served (FCFS) service
discipline and general distribution function of service time at the servers,
has proved to be resilient to all approaches except for approximate solution
techniques. In this chapter our attention is limited to this type of queueing
network.
Several approximation methods have been suggested for its treatment.
On the one hand there are diusion approximations [9, 10, 20] applicable
to two-station networks or to general queueing networks [11, 12, 13, 16, 22]
and on the other hand we have iterative techniques [5, 23]. The convergence
of the latter to the exact solution is not an established fact and we know
that the former tend, in certain simple cases, to the exact solution.
Most of the work published in the literature has concentrated on
evaluating the joint probability distribution of queue lengths for all the
queues in a network, but it is seldom possible to make use of this complete
information. In measurements on computer systems it is dicult enough to
collect data on the performance of a single resource, and the measurement
of joint data for several resources could become very time- and space-
consuming. The same can be said of simulation experiments where the task
of computing condence intervals for estimated joint statistics becomes
impractical. Furthermore, when it comes to computing average response
times or queue lengths it suces to know the average response time
encountered in each individual queue. Therefore it would suce in many
cases to be able to compute with satisfactory accuracy the probability
distribution for the queue length at each individual resource.
165
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
There are two ways in which we may intuitively understand the basis
for diusion approximations. The rst uses a numerical analysis analogy,
while the second calls upon the central limit theorem.
2
f (x, t) b f (x, t) + f (x, t) = 0 (6.1)
t x 2 x2
where f (x, t) is a function of space, the x variable, and of time t. f (x, t) is
chosen to be, for each t, the probability density function of a non-negative
random variable X(t):
Equation (6.1) can be solved if an initial condition (f (x, 0) for all values
of x 0) and a boundary condition (conditions which must be satised
by f (0, t) and f (0, t)/t) are provided. We shall consider the following
boundary condition given in terms of P (t) a probability mass, function of
time, located at the boundary point x = 0:
d 1 f (x, t)
P (t) = cP (t) + lim bf (x, t) +
dt x0+ 2 x (6.2)
f (0, t) = 0.
Equations (6.1) and (6.2) are to be viewed, for the moment, simply as
formal relations. The interpretation in terms of queueing phenomena can
be obtained either in terms of their discretisation (as will be done here), or
via the central limit theorem as in section 2.2.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
so that
M
M
pi (t) f (x, t)dx + P (t).
i=0 0
c = /2, b = 0, = 1
we see that (6.5) are the ChapmanKolmogorov equations for the M/M/1
queue with arrival rate equal to the service rate (see Chapter 1). Thus
we can conclude that for this special case ( = ), the diusion equation is
approximated by the equation for the M/M/1 queue and vice versa.
If we seek the stationary solution of (6.1), (6.2) or (6.5), then these
equations are approximations of each other under general conditions.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
= /2, = (/2) b, or b =
c = (1 + b/) = 2 /.
where P and f (x) denote the stationary solution. Similarly, for (6.5) we
have
pi = (/)i (1 /)
P p0 ; f (i) pi /, i 1.
or
E[Q(t)] = ( )T = bT
and
2
f (x, t) b f (x, t) + f (x, t) = 0
t x 2 x2
where {X(t), t 0} is the continuous-path stochastic process approximat-
ing the number in queue.
Since the approach was initially intended for heavy trac conditions
it is also assumed that the lower boundary at x = 0 for the process
{X(t), t 0} should act as a reecting boundary. This last assumption
implies that no probability mass can collect at x = 0.
From (6.1) we may write (for f = f (x, t))
f 2f f
dx = b f + dx = bf + .
0+ t 0+ x 2 x2 2 x 0+
The left-hand side must be zero because the total probability mass is one,
and no probability mass collects at x = 0. Therefore, for all t 0,
f +
bf (0+ , t) = (0 , t)
2 x
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
f (x) = ex , x 0
= 0, x<0
p(0) = 1
i1 ,
p(i) = (1 ) i1
for = e .
where
1 if i = 1
ai =
b1 . . . bi1 if i > 1, 0 < bi 1.
where
1 if i = 1
Ai =
B1 . . . Bi1 if i > 1, 0 < Bi 1.
At the end of the holding time at the upper boundary the particle
jumps back instantaneously to a random point in ]0, M [ whose position
is determined by the probability density function f2 (x).f1 (x) and f2 (x)
may be taken to be functions of the instant at which the jumps occur.
Notice that
n
ai
E[h] = ;
i=1 i
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
similarly,
m
Ai
E[H] = .
i=1
i
= (E[h])1 , = (E[H])1 .
1 2
Ax,t f = f bf + f
t x 2 x2
1
Cx,t f = bf + f.
2 x
Also, let Pi (t), 1 i n, be the probability that the particle is in the
i-th stage of the holding time at the lower boundary at time t while Qi (t),
1 i m, is the probability that it is in the i-th stage of the holding time
at the upper boundary at time t. The equations describing the evolution of
the particle are
n
m
Ax,t f + i (1 bi )Pi (t)f1 (x) + i (1 Bi )Qi (t)f2 (x) = 0 (6.6)
i=1 i=1
d 1 P1 (t) + C0,t f if i = 1
Pi (t) = (6.7)
dt i Pi (t) + i1 bi1 Pi1 (t) if 1 < i n
d 1 Q1 (t) CM,t f if i = 1
Qi (t) = (6.8)
dt Q (t) + B Q (t) if 1 < i m
i i i1 i1 i1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
where
1
C0,t f = lim bf + f
x0 2 x
1
CM,t f = lim bf + f
xM 2 x
= (E[h])1 , = (E[H])1 .
Dene P (t) as the probability that the particle is at the lower boundary
at time t, and let Q(t) be the corresponding quantity for the upper
boundary:
n
m
P (t) = Pi (t), Q(t) = Qi (t).
i=1 i=1
m
d
Q(t) = i (1 Bi )Qi (t) CM,t f. (6.10)
dt i=1
m
+ i (1 Bi )Qi (t) f2 (x)dx (6.11)
i=1
which states that the rate of change of the probability mass in is equal
to the rate of ow of the probability mass out of (the rst term on the
right-hand side of (6.11)) plus the rate of ow into from x = 0 and from
x = M (the second and third terms, respectively, on the right-hand side). In
order to deduce (6.7), notice that for 1 < i n we may write for any t 0,
since the time the particle spends in any one of the stages of the Cox
distribution is exponentially distributed; by collecting terms, dividing both
sides by t and taking t 0, this yields (6.7) for 1 < i n in the
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
P1 = 1
1 C0 f, Pi = (i1 bi1 /i )Pi1 , 1<in
so that
ai
Pi = i b1 . . . bi1 C0 f = C0 f, 1 < i n.
i
Therefore,
n
P = Pi = 1 C0 f ;
1
Q = 1 CM f.
But
n
n
i (1 bi )Pi = ai (1 bi )C0 f = C0 f
i=1 i=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
and similarly
m
i (1 Bi )Qi = CM f.
i=1
Since these equations depend on E[h] and E[H] only we have proved
that the stationary probabilities f (x), P, Q are independent of the higher
moments of h and H.
This result shows that if we are interested in approximating the
stationary queue length probability distribution using the instantaneous
return process, it suces to use a model where h and H are exponential
This is the assumption we will make in the sequel.
b =
= 3 Va + 3 Vs = Ka2 + Ks2 .
queue in a busy period to the rst arrival of the next busy period. If the
arrival process is Poisson it is natural to take = (E[h])1 . However, if the
arrival process is not Poisson then the interarrival time distribution and
the distribution of h need not be the same. Let us consider the case where
(E[h])1 = = .
The instantaneous return process approximation in stationary state for
the GI/G/1 queue is represented by the equations obtained from (6.12),
(6.13):
f 1 2t
b + 2 + P (x 1) = 0 (6.15)
x 2 x
P = C0 f. (6.16)
Notice that the term f1 (x) in (6.12) has been replaced by the Dirac density
function concentrated at x = 1, (x 1). This represents the fact that when
an arrival occurs, the queue length jumps instantaneously from x = 0 to
x = 1. In addition to (6.15) and (6.16) we also use f (0) = 0 and
P+ f (x)dx = 1.
0+
R = /( + ). (6.19)
L = 1 (1 Ks2 )
L
2
L)/L
so that the relative error (L tends to zero as 1.
Fig. 6.1.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
where = /, and
K = (1 2 e(M1) )1 .
Fig. 6.2. Maximum percentage relative error of diusion approximation for closed two-
server system. , exponential IOD; , constant IOD; , condence interval for
exponential IOD; , condence interval for constant IOD.
P = 1 .
which is
i = [1 ]
1]
p1 (i) = [ i1
p1 (0) = 1
where = e ; notice that this is identical to (6.5). The average queue length
obtained is then
L1 = ip1 (i) = /(1 )
i=1
Property (ii) is important since it states that as Ks2 increases, the relative
error depends only on ; but the factor 2(1 ) will be unacceptably
high for small values of . Property (i) is a general property of diusion
approximations: the relative error tends to zero under heavy trac
conditions.
Let us examine how L1 behaves when Ks2 is large and small, i.e. when
L2 is a poor approximation. With Ka2 = 1, we have
which is
= exp[(2/Ks2 ) (1 /Ks2 )(1 )]
for Ks2 1 , or
= 1 2(1 )/Ks2
so that
L1
= K 2.
2(1 ) s
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
( 2)
= Ks2 + .
2 2(1 )
Thus, for
1,
just as for L2 .
Consider now the case where = (1 )
1. This analysis has been
carried out in [25]. We will have
2
2 1 2
= 1 + + O(3 )
+ Ks2 2 + Ks2
so that
( + Ks2 ) (1 ) 2
L1 = 1+ + O( )
2(1 ) + Ks2
and
2 1 2
L1 LPK = Ks + + O( ) .
2 + Ks2
Therefore:
(iii) lim (L1 LPK )/LPK = 0, and for = (1 )
1, we have
1
1
(iv) lim (L1 LPK )/LPK = + O(2 )
Ks2
while for L2 , we have (ii); thus we see that for close to 1, L1 has a relative
accuracy which is twice as good as that of L2 .
In fact, the form of the PollaczekKhintchine formula suggests a new
approximation which was noticed in [10] and further developed in [6].
Instead of choosing as has been done above, suppose that we take
2(1 )
= .
(Ka2 + Ks2 )
We may then derive
1 (Ka2 + Ks2 )
L2 = 1 = 1 +
2(1 )
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
using the same form for the probabilities as given in (6.25) except that we
replace by and by :
= e .
Table 6.1. Exact and approximate average queue lengths of the M/G/1
queue for = 0.8
Table 6.2. Approximate average queue length for the E2 /H2 /1 system
compared with simulation results (95% condence intervals) for Ka2 = 0.5
The results of Table 6.3 use the well-known fact that in general the
stationary solution of the GI/M/1 queue is given by (see, for instance, [14]):
p0 = 1
pi = (1 )i1 , i1
= A ( )
and A (s) = 0 esx dA(x), where A(x) is the interarrival time distribu-
tion. Therefore, the average stationary queue length of the G1/M/1 queue
is /(1 ).
The various analytical results and numerical examples are evidence
to the eect that the heuristic modication should be chosen. The
corresponding diusion parameters are b = and = Ka2 + Ks2 .
These will be the values retained in the following sections, so that the
discretised approximation which we shall use is
p0 = 1
1
p1 = 1 (
1) ,
pi = (1 )2 i1 , i 2 (6.26b)
with = e .
Fig. 6.3. General open queueing network with rst-in-rst-out service discipline.
which is unique under these assumptions (see Chapter 3). Then ei is the
expected number of visits which a customer of the network will make to
station i. The arrival rate of customers to station i is i = 0 ei at steady-
state; also the steady-state probability i that station i contains at least
one customer is given by
0 ei
i = if ei < i
i
where
1
i = t dFi (t)
0
is the average service time for a customer at station i. This fact can be
easily established rigorously; one way is to treat an open network of this
kind as a limiting case of a closed network when one station is saturated
and to apply the work-rate theorem (see Chapter 3).
The approach we develop in this section is based on the following
assumption, which in general is unjustied: the departure process from
any station in the open network is a renewal process, i.e. times between
successive departures are independent and identically distributed. We shall
make use of this assumption in order to compute the rst two moments of
the interdeparture time distribution, although it is in general not satised.
This assumption is valid in the open network with Poisson arrivals and
exponentially distributed service times. It is also valid for the output of
station i when 0 ei /i 1, or when all j = 0. Let Ci , 1 i n, be
the squared coecient of variation of the interdepartures times at station
i, and denote by Ai the interarrival time, Si the service time, Ai the idle
time, and by i the interdeparture time. We shall dene C0 = K02 in order
to maintain a uniform presentation.
For t large enough, and assuming that the output processes from each
individual queue are independent, the total number of arrivals to station i
in the interval [0, t] will be normally distributed with mean i t and variance
n
[(Cj 1)pji + 1]j pji t.
j=0
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
Here we have used the fact that the sum of independent normal random
variables is normal with variance being the sum of individual variances. In
the usual diusion equations for approximating the length of each queue
(6.15), (6.16), the following parameters will be chosen, following (6.26) (see
section 6.2.6):
bi = i i , i = i , i = i /i
n
(6.27)
2
i = i i Ki + [(Cj 1)pji + 1]j pji
j=0
where the subscript i refers to the parameters of the equations of the i-th
queue, and Ki2 is the squared coecient of variation of service time at the
i-th queue.
In order to complete the development we must obtain Ci , 1 i n. We
shall assume that i is a service time Si with probability i or an interarrival
time plus a service time Ai + Si with probability (1 i ). We then have
E[i ] = i 1 1 1 1
i + (1 i )(i + i ) = i
E[i2 ] = (1 2 2 2
i ) (1 + Ci ) = E[Si ] + (1 i )(E[Ai ] + 2E[Ai ]E[Si ])
so that
Finally, it is the
n
[(Cj 1)pji + 1]j pji = 3i (E[A2i ] (1 2
i ) ).
j=0
fi (xi , t) is the density function approximating the length of the i-th queue
and Pi (t) is the probability that the i-th queue is empty. The stationary
solution will be obtained as
i (e i 1)ei xi , xi 1
fi (xi ) =
i (1 ei xi ), 0 xi 1
Pi = 1 i
Li = i [1 i /2bi ]
The system of equations (6.28) is then solved with the modied values pij
and K 2 . Notice that the arrival rate to a queue is modied only if pii = 0
i
in the original network; however, the value of the load factor i = i /i
is unchanged since i becomes i (1 pii ) and i becomes i (1 pii ).
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
The queue length is also preserved, since service times are being replaced
by longer service times (see (ii)) which are the sum of a geometrically
distributed number of service times corresponding to the feedback of
customers to the queue. That is, Si is being replaced by
l
Si = Sik
k=1
Fig. 6.4.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
being [(1 21 )/2 + K12 /21 ]; the value of C1 /2 obtained from equation
(6.28) is exactly this value.
Example 6.1
In order to illustrate the degree of accuracy which can be obtained from the
approximation techniques for open single-customer class queueing networks
which we have presented in this section, we shall apply the preceding
results to the system model presented on Fig. 6.5. The model which is
shown was introduced in [1] in order to evaluate the performance of an
interactive system. In this model, jobs arrive at the system in a Poisson
stream of rate 0 ; after passing through server 1 (which represents a
central processing unit) they either leave the system with probability
(1 1 ) or enter the queue of server 2 (representing an input-output
device). A customer will either enter once again the queue of server 2, with
probability 2 , or proceed to server 1 after nishing its service at server 2.
The following simulations and numerical results have been obtained by
Dinh [7]. The results shown in Table 6.5 provide a comparison, with
respect to simulation results, of the accuracy of the method developed
earlier in this section as well as of the approach of Reiser and Kobayashi
which we have summarised in section 6.3.2. Condence intervals for the
values estimated from the simulation experiments have not been provided.
However, the precise conditions under which the simulations were carried
out are shown in Table 6.4 for the ve simulation runs. These data indicate
that it was impossible in any of the simulations to obtain the same arrival
spi-b749
Experiment Queue
No. 1 2 0 K0 No. 1
i Ki Simulated Computed Simulated Computed
1 0.510 0.503 0.512 0.941 1 0.91123 0.427 1.0489 1.045 0.957 0.953
2 0.84000 0 1.078 1.073 0.905 0.901
2 0.509 0.499 0.410 0.944 1 0.91591 0.423 0.835 0.836 0.764 0.766
2 0.84000 0 0.848 0.849 0.712 0.713
3 0.516 0.506 0.342 0.945 1 0.91443 0.414 0.707 0.706 0.646 0.646
2 0.84000 0 0.738 0.738 0.620 0.620
9in x 6in
4 0.512 0.502 0.293 0.967 1 0.90436 0.432 0.601 0.602 0.544 0.544
2 0.84000 0 0.619 0.619 0.520 0.520
5 0.504 0.507 0.257 0.952 1 0.91094 0.422 0.519 0.519 0.472 0.473
2 0.84000 0 0.530 0.531 0.444 0.446
b749-ch06
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
spi-b749
Experiment Queue G/P G/P
No. No. Simulation G/P Simulation G/P R/K Simulation R/K No mod. with mod.
9in x 6in
4 1 1.070 0.785 1.189 0.898 0.729 1.046 0.815 0.821 0.818
2 2.858 0.663 3.008 0.860 0.603 1.005 0.545 0.502 0.781
5 1 1.184 0.824 1.284 0.912 0.728 0.758 0.632 0.619 0.617
2 3.341 0.740 3.438 0.890 0.599 0.731 0.453 0.382 0.589
b749-ch06
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
1 = rij
= (c1 , d1 )
k+1 = (ck+1 , dk+1 ) = rd k ,j , if dk = j
L = k, if dk = m.
i dij
=
or the proportion of packets which enter link l after having entered link k.
The k and pkl obtained from (6.30) and (6.31) can now be used in (6.28)
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
where k = k /k , 1
k is the average transit time of a packet along link k.
If the usual formula [15] using Jacksons theorem had been used, we would
have had (see Chapter 3):
Example 6.2
Let us now apply these results to a numerical example given in [4]. The
network is shown in Fig. 6.6 and the external arrival process is Poisson
with rates:
1 = 6;
2 = 8.25;
3 = 7.5;
4 = 6.75;
5 = 1.5.
Packet lengths are assumed to be constants so that Kk2 = 0 for all links
1 k 12, and the packet length is 1000 bits. Links 1, 2, 7, 8, 11, 12 have
a data-transmission capacity of 4800 bits/second, while links 3, 4, 5, 6, 9,
10 have a capacity of 48,000 bits/second. Therefore
1 = 2 = 7 = 8 = 11 = 12 = 4.8
3 = 4 = 5 = 6 = 9 = 10 = 48.
0 3 3 3 2
4 0 5 5 4
R=
6 6 0 9 8 .
10 10 10 0 12
1 1 7 11 0
The system as described was simulated until 6000 packets were received
at their destinations. Although the simulation results have not been
analysed in order to compute statistical condence intervals, this simulation
experiment is comparable in duration to a measurement session on a real
computer network, the main point being that the diusion approximation
is capable of making predictions which are as accurate as simulation of such
systems. The results obtained are given in Table 6.6.
We see that for buer queues which are lightly loaded, Jacksons
formula (6.33) yields results which are of the same degree of accuracy as
the diusion approximation. However for link 2, which is heavily loaded,
the diusion approximation is considerably more accurate.
Let us complete this section by deriving another formula which is of
interest in the analysis of packet-switching networks. A useful performance
measure is the average source-destination transit delay T(i,j) for the source-
destination pair (i, j). For xed routing this corresponds simply to the
Table 6.6. Average buer queue lengths for the packet-switching network
of Fig. 6.6
average time to traverse the path (i, j) which corresponds to (l, m).
Therefore
T(i,j) = Lk /k . (6.34)
k(i,j)
The important case of networks with nite storage capacity will not be
analysed here, and the interested reader is referred to [3] for results on this
subject.
i = ei
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
where
2(i i )
i = . (6.37)
i K (i) + i Ki
The approximation proposed for the joint probability distribution is
n
p(m1 , . . . , mn ) = pi (mi ) (6.38)
i=1
for the n queues in an open network. Obviously the result holds only if
i i , 1 i n, which is the usual stability condition.
For a closed network, the following treatment has been suggested.
Suppose i is the utilisation of server i; the joint probability distribution is
taken to be (for M customers in the network):
n
p(m1 , . . . , mn ) = G pi (mi ) (6.39)
i=1
where
1 i , mi = 0
pi (mi ) = (6.40)
i (1 m
i )i
i 1
, mi = 1, 2, . . . , M
i = Xi K /Xk i (6.41)
(i) a stream of arrivals to the network which is a renewal process: its rate
is 0,r and the squared coecient of variation of interarrival times is
K0,r ;
(ii) a general service time distribution function Fri (x) at the i-th service
station, 1 i n. 1 i,r will be its average and Ki,r its squared
coecient of variation.
Let us denote i,r = i,r /i,r : it can be viewed as the load imposed by class
r customers on station i. We shall dene
R
i = i,r (6.43)
r=1
and
R
i = i,r , i,r = i,r /i for 1 i n, 1 r R,
r=1
R
(6.44)
0 = 0,r ,
r=1
We also obtain:
Ci + 1 = 2i E[Si2 ] + (1 i )(Gi + 1) + 2i (1 i ).
The service time Si is Si,r if it is the service of a class r customer (i.e. with
probability i,r = i,r /i ). Therefore
R
E[Si2 ] = 2
E[Si,r ]i,r
r=1
R
= 2
(Ki,r + 1)(1 2
i,r ) i,r /i
r=1
so that
R
Ci + 1 = i i,r 1 2
i,r (Ki,r + 1) + (1 i )(Gi + 1 + 2i ). (6.46)
i=r
Using an argument similar to the one used in section 6.3, assuming that
the output processes of the n queues are mutually independent renewal
processes, we can write that the variance of the number of arrivals at the
i-th queue in a long interval (0, t) will be
n
3i {E[A2i ] (E[Ai ])2 }t = i Gi t
= [(Cj 1)pji + 1]j pji t
j=0
where pji is the probability that a job leaving station j enters station i, or
R
R
pji = j,r pj,r;i,r (6.48)
r=1 r =1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
Using (6.46) and (6.47) we now obtain the system of n linear equations
for the Ci , 1 i n:
R
n n
j j 2
Ci = i i,r 1 2
i,r (Ki,r +1)+ (1pji ) pji 22i + p Cj . (6.50)
i=1 j=0
i ji
j=0 i
Begin
Step 1 Obtain the i,r , 1 i n, 1 r R from the linear system of nR
equations (6.42).
Step 2 Compute i , i , i,r from (6.43), (6.44).
Step 3 Use (6.48) to obtain the pji , 0 j n, 1 i n.
Step 4 Obtain the Ci , 1 i n, by solving the n linear equations (6.50)
using C0 from (6.49).
Step 5 Compute Gi , and G1,r , 1 i n, 1 r R, from (6.47) using the
result of Step 4, and using (6.51).
End.
fi fi 1 2 fi
bi + i 2 + i Pi (t)(xi 1) = 0
t xi 2 xi
d 1 fi
Pi (t) = i Pi (t) + lim+ bi fi + i
dt xi 0 2 xi
with lim fi (xi , t) = 0, and where Pi (t) is the probability that the queue
xi 0+
length is xi = 0 at time t. The parameters for the diusion process are
chosen to be
bi = i i , i = i i Gi + 3i Vi (6.52)
Pi = 1 i /i = 1 i (6.54)
i [1 ei xi ], 0 xi 1
fi (x) = (6.55)
i [1 ei ]ei xi , xi 1
From the distribution for the total number in queue we will now work back
to the distribution of the number of customers of each class in queue. We
proceed as follows. Discretise the probability density function fi (x) by using
(6.26b). pi (ni ) will be the discrete approximation to the stationary queue
length distribution at station i. Let pi,r (li ) be the probability of nding
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
Wi,r = Li /i 1
i
and
Ti,r = Li /i 1 1
i + i,r . (6.58)
Therefore
n
Tr = Ti,r i,r /0,r (6.59)
i=1
since a class r customer will visit station i on average i,r /0,r times.
A detailed validation of this models predictions is given in [7] for a
computer system with two job classes. The accuracy seems to be very good.
Comparisons with simulation results reported in [7] yield a relative error of
less than 10% in average queue lengths for each class.
6.5. Conclusion
Many important practical cases of large-scale computer systems are too
complex to be represented exactly by a mathematical model. Even when a
precise mathematical model can be constructed, the analyst is faced with
a problem of dimension. Models with a number of states proportional to
106 are easy to obtain, but program packages capable of solving Markov
chains of this dimension are not yet available. Often the mathematical
models which arise from computer systems have properties which make
them particularly dicult to handle numerically; for instance, the time
constants related to various parts of the system will vary widely leading to
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
References
1. Anderson, H. A. and Sargent, R. (1972). The statistical evaluation of the
performance of an experimental APL/360 system. In Statistical Computer
Performance Evaluation (W. Freiberger, Ed.), pp. 7398. Academic Press,
London.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch06
Chapter 7
Approximate Decomposition and Iterative
Techniques for Closed Model Solution
7.1. Introduction
In general, the computer system analyst has to adapt the panoply of tools at
his disposal to the specic problem at hand, and often his problem will not
t exactly into any available framework. If the analyst has a mathematical
orientation, and enough time available, he may attempt an original solution
method. If he is pressed for time, or if he is not mathematically inclined,
he will tend to program a simulation of the system he has to analyse if the
programming and computer time can be aorded.
There is, however, a third approach he might take: the use of numerical
approximation which in some cases provides only a rst-order approxima-
tion, and which in others provides highly accurate results. The diusion
approximation developed in Chapter 4 is an example of this approach. In
this chapter we will examine a set of approximations which retain, contrary
to diusion approximations, the discrete nature of the model. They will all
be based on a similar approach to a heuristic iterative solution of the steady-
state Birth and Death equations. However, in certain cases, a formal
justication will be available on the basis of problem structure while in
other cases the only justication will be the intuitive appeal of the approach
and its similarity with techniques used in other areas of applied science.
211
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch07
under the eect of the rest of the system viewed as an aggregate. In certain
cases, we consider several subsystems interacting with each other and with
the rest of the model. The simplest such structure is shown in Fig. 7.1,
where the subsystem is a single queue while the aggregate contains the
remaining queues of the system. An iterative solution technique based on
an aggregate/subsystem decomposition will often iterate between several
dierent decompositions of the type shown in Fig. 7.1 in order to reach an
approximate solution. The stationary solution will then be framed either in
terms of marginal distributions for each subsystem or as a product of the
marginal distributions when the stationary distribution for global network
state is desired.
In order to illustrate this common heuristic approach, we will rst apply
it to a class of queueing networks for which it yields an exact result: a closed
network of exponential queues with FIFO service discipline and a single
class of customers.
n = (n1 , . . . , nN )
Finally, we relate the input ow to any queue to the output ows from all
the queues:
N
i = j Pji
j=1
Now, assuming that the aggregate and the subsystem (the N -th queue)
interact very infrequently, we can write using (3.37)
ei GN 1 (K nN 1)
p(n1 , . . . , nN 1 , nN ) =
n >0
i GN 1 (K nN )
PN 1 i
nj =KnN
or using (7.1):
pN (nN ) eN GN 1 (K nN )
= . (7.3)
pN (nN 1) N GN 1 (K nN + 1)
This procedure yields in fact the exact result for a Jackson network.
Using (3.37) and the argument leading to (3.38) we can see that
(eN /N )nN GN (K nN ) (eN /N )nN +1 GN (K nN 1)
pN (nN ) = .
GN (K)
But using (3.35) we obtain
nN
eN GN 1 (K nN )
pN (nN ) =
N GN (K)
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch07
which obviously satises (7.3). Therefore (7.2) is in fact the exact form of
the solution.
with (7.4) instead. We maintain the usual assumption that the system
is strongly connected (irreducible and aperiodic), i.e. each state can be
reached from any other state with non-zero probability in a nite amount
of time.
Let be a partition of the set of service centres {1, . . . , N }: =
{ 1 , . . . , l }. We shall say that two state vectors n and n are -equivalent
if and only if
(j) ni = ni .
i j i j
= { 1 , . . . , k }.
where n -/ -n means that n and n are not -equivalent. Inequality (7.5)
stipulates that the transitions between states which are -equivalent are far
more likely than others, and hence far more frequent. Thus, if a queueing
network is NCD on a partition of its service centres, changes in the
number of customers in each element of will be relatively infrequent with
respect to state transitions which do not modify that number.
We shall say that a partition = { 1 , . . . , l } is non-trivial if
1 < l < N . An element of a non-trivial partition will be non-trivial if it
contains more than one service centre. Henceforth we will consider a non-
trivial partition; let k be one of its non-trivial elements, and each element
i of corresponds to a set of -equivalent states:
Q may be written as
Q = D + E (7.6)
Since both Q and D are stochastic matrices, it follows that row sums of
E will be zero. Furthermore, |e(n, n )| 1.
It is clear from (7.7) that by a permutation of rows and columns, D
may be rewritten in block diagonal form as:
(7.9)
q = qQ (7.10)
with
q(n) = 1
n
d = dD (7.11)
will be used to approximate q. In fact (7.11) does not have a unique solution,
even when the condition
d(n) = 1
n
q = qQ = qD + qE (7.12)
q=d+ (7.13)
d + = (d + )D + qE
or
(I D) = qE. (7.14)
Example 7.1
Let
1 3 1
8 8 2
1 1 1
A=
4
.
4 2
1 1 1
3 3 3
A is lumpable on the partition = {(1, 2), (3)} and
1 1
2 2
A =
2 1.
3 3
Example 7.2
The matrix D dened by (7.4) and shown in (7.9) is lumpable on
= {1 , . . . , k } and
1 0
1
D = . .
..
0 1
is the k k unit matrix.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch07
d = dD.
= [Q E] + qE. (7.17)
a1 = a1 Q = qE
(7.18)
ai+1 = ai+1 Q ai E, i 1.
a1 .
q(n, n ), otherwise.
and by subtracting the sum of the quantities added from the diagonal
elements; is taken to be
k
= max ij . (7.23)
i
j=1
j=i
This guarantees that all elements of U will be less than one in absolute
value. The following points should be noticed:
be the vector of stationary probabilities associated with Q:
Let q
=q
q
Q.
Write
+
q=q
+=q
q + Q
+ qU
or
+ qU = [Q U] + qU.
= Q (7.24)
We write
= bi i (7.25)
i=1
= b1 (7.27)
q
= d + b1 + a1 (7.28)
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch07
if we write
q d) =
( ai i
i=1
K1
2 (K i + 1)
p1 (K1 ) = p2 (K K1 ) = p1 (0)
i=1
1 (i)
where
K1
K
1
2 (K i)
p1 (0) = 1 + .
i=1
1 (i)
K1 =1
N
ej pji = ei
1
N
Uj,m pji = Ui,m
j=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch07
N
(i) (1 )K i,m (1 + )K,
Q and
i=1
N N
1 1
(ii) (1 ) j,m i,m (1 + ) j,m ,
N j=1 N j=1
Fig. 7.4.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch07
If
then the procedure stops at the m-th iteration: the quantities computed
in step (1.2) are considered to be a satisfactory approximation to R.
Otherwise start at step (1) with m replaced by m + 1 and i,m+1 computed
in step (1.5).
Step (2): If (i) is not satised but (ii) is satised, go to step (2.2);
otherwise proceed to step (2.1).
Step (2.1): Compute, for 1 i N ,
N
i,m+1 = i,m i,m N j,m .
j=1
If
The role of the steps of this iterative technique merits some explanation.
All parts of steps (1) are used to compute an approximation from Rm
to the related quantities of the original network R. Step (2) constructs
modications to Rm which will be incorporated in Rm+1 .
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch07
References
1. Chandy, K. M., Herzog, U. and Woo, L. (1975). Parametric analysis of general
queueing networks. IBM Res. and Dev., 19, 3642.
2. Chandy, K. M., Herzog, U. and Woo, L. (1975). Approximate analysis of
general queueing networks. IBM J. Res. and Dev., 19, 4349.
3. Courtois, P. J. (1972). On the Near-Complete Decomposability of Networks of
Queues and of Stochastic Models of Multiprogramming Computer Systems.
Computer Science Report, CMU-CS-72, III, Carnegie-Mellon University,
Pittsburgh, Pennsylvania.
4. Courtois, P. J. (1977). Decomposability: Queueing and Computer System
Applications. Academic Press, New York.
5. Marie, R. (1978). Methodes iteratives de Resolution de Reseaux de files
dAttente, Doctoral Thesis, Universite de Rennes.
6. Shum, A. V. and Buzen, J. P. (1977). A method for obtaining approximate
solutions to closed queueing networks with general service times. In
Modelling and Performance Evaluation of Computer Systems (H. Beilner
and E. Gelenbe, Eds), pp. 201220. North-Holland, Amsterdam.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch07
Chapter 8
Synthesis Problems in Single-Resource
Systems: Characterisation and Control
of Achievable Performance
231
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
W = (W1 , W2 , . . . , WR )
S W.
Fig. 8.1.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
Since VS (t), and hence E[VS (t)], is independent of S for every t, we have
the following basic result.
Theorem 8.1 (General Conservation Law). For any single-server
queueing system in equilibrium there exists a constant V, determined only
by the parameters of the arrival and required service times processes,
such that
VS = V (8.1)
where VS (r) is the expected steady-state virtual load due to jobs of class r
(the sum of the average remaining service times of all class r jobs in the
system at a random point in the steady-state) under scheduling strategy S.
For a given r, the value of VS (r) depends on S in general (e.g. if the priority
of class r jobs is increased, VS (r) is likely to decrease). Theorem 8.1 asserts
that the vector (VS (1), VS (2), . . . , VS (R)) always varies with S in such a
way that the sum of its elements remains constant. Note that the truth of
this statement does not rely on any assumptions about interarrival times,
service times or independency between them.
Intuitively, the average virtual load due to class r is related to the
average number of class r jobs in the system, and hence to the average
response time for class r jobs. The general conservation law (8.2) should
therefore imply a relation among the elements of the response time vector
W. However, in order to render such a relation explicit we have to make
more restrictive assumptions regarding the nature of the demand processes
and the complexity of the scheduling strategies.
A scheduling strategy is a procedure for deciding which job, if any,
should be in service at any moment of time. It takes as input any
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
information that is available (the time of day, the types of job in the system,
their arrival instants, the amounts of service they have received, etc.) and
returns the identier of the job to be served, or zero if the server is to
be idle. The only restriction imposed so far has been that the procedure
returns zero if, and only if, there are no jobs in the system. Now we add
two more restrictions:
(i) every time the server becomes idle the procedures memory is cleared;
the scheduling decisions made during one busy period are not based on
information about previous busy periods;
(ii) only information about the current state and the past of the queueing
process is used in making scheduling decisions; thus, it is possible to
discriminate among jobs on the basis of their expected remaining service
times (since their types and attained service are known), but not on
the basis of exact remaining service times.
VS (r) = Nr /r , r = 1, 2, . . . , R,
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
VS (r) = r Wr , r = 1, 2, . . . , R (8.3)
Theorem 8.2. When the required service times are distributed exponen-
tially, there exists a constant V determined only by the interarrival time
distributions and by the parameters r , such that
R
r Wr = V, (8.4)
r=1
job in service, given that it is of class r, is equal to the average residual life
r of the class r service time. r is given by
1
r = M2r r , r = 1, 2, . . . , R, (8.5)
2
where M2r is the second moment of the class r service time distribution
(see (1.66), Chapter 1).
We can now write, for the average virtual load due to class r jobs,
r = 1, 2, . . . , R,
1 1
VS (r) = nr + r r = r wr + r r = r Wr r r . (8.6)
r r
= 1 + 2 + + R ,
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
where
R
w0 = r r (8.9)
r=1
R
R
1 r
r Wr = . (8.10)
r=1
1 r=1 r
R
R R
r
w0 1 w0
r Wr = + r r = + (8.11)
r=1
1 r=1 r 1 r=1 r
where
Does such a minimising strategy S exist? If pre-emptions are allowed,
the answer is clearly yes: any strategy which gives pre-emptive priority to
g-jobs over non-g-jobs can be taken as S , since all these strategies eliminate
the g-intervals completely. (8.12) can then be rewritten as
VS (r) V g , for all S (8.13)
rg
Note that the Poisson input assumptions were used only in order to
write a closed-form expression for the right-hand side of (8.13); if we leave
it as V g , Theorem 8.4 will continue to hold.
The situation is less straightforward if one is restricted to non-pre-
emptive scheduling strategies only. Now, if g is a proper and non-empty
subset of {1, 2, . . . , R}, the inuence of the jobs whose classes are in
{1, 2, . . . , R}-g cannot be eliminated completely. There is no scheduling
strategy which minimises the g-intervals for every realisation of the demand
processes. However, for a given realisation, the strategy which minimises
the g-intervals has to be one that gives head-of-the-line priority to g-jobs
(eliminating all g-intervals except, perhaps, those at the start of g-jobs
busy periods). Therefore, only such a priority strategy can minimise the
steady-state average virtual load due to g-jobs, VSg . Making the appropriate
assumptions and using (8.6) we can rephrase the above statement thus: in
order to minimise the linear combination rg r Wr it is necessary to give
non-pre-emptive priority to g-jobs.
Now suppose that the input streams are Poisson. If the g-jobs have
non-pre-emptive priority, the only way their average response time can be
inuenced by the non-g-job is through the probability that an incoming
g-jobs nds a non-g-jobs in service. But with Poisson arrivals, that
probability is independent of the scheduling strategy (see Chapter 1).
Hence, if g-jobs have non-pre-emptive priority, VSg is independent of the
order of service among the non-g-jobs. It is also independent of the order
of service among the g-jobs because, once the g-jobs have started being
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
served, there are no g-intervals until the end of the busy period. Thus, the
minimal value of rg r Wr can be obtained by lumping all g-jobs in one
class, all non-g-jobs in another class, and giving head-of-the-line priority
to the g-jobs. Performing these calculations yields a generalisation of
conservation law 2:
1
r Wr w0 r 1 r + (r /r ) (8.15)
rg rg rg rg
In the next section, the relations derived here will lead to a charac-
terisation of the sets of achievable performance vectors. Before proceeding,
however, we should take another look at the assumptions that have been
made and at the possibilities for relaxing them.
First, we shall consider the scheduling strategies. It is evident that
if the strategies are not required to be work-conserving, Theorem 8.1 and
all that follows from it will hold no more. The necessity of condition (ii) for
the special conservation laws is less obvious (that condition is very rarely
mentioned in the literature) but it, too, turns out to be unavoidable. We
shall give examples of both pre-emptive and non-pre-emptive scheduling
strategies where exact service times are known in advance and where (8.10)
and (8.11) do not hold.
Could we drop the exponential service times assumption and still allow
pre-emptions? The answer is again, alas, no. Neither (8.10) nor (8.11) are
satised in the case of the classic pre-emptive priority disciplines with
general service times.
Finally, we know that the Poisson inputs assumption is not necessary
for the validity of Theorem 8.4. What is not known is whether Theorem 8.5
continues to hold (perhaps with dierent constants on the right-hand side
of the inequalities) if that assumption is relaxed.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
H1 H1 . (8.16)
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
Fig. 8.2.
H2 H2 . (8.17)
W1 = W1 (i1 , i2 , . . . , iR ).
W2 = W2 (i1 , i2 , . . . , iR ).
where one of the gj is the set {1, 2, . . . , R} and the others are proper, non-
empty and dierent subsets. Using the notation ar = r /r , ag = rg ar
and g = rg r , we rewrite these equations as
r Wr = agj /(1 gj ); j = 1, 2, . . . , R. (8.18)
rgj
where gjk = gj gk ; the last term is zero by denition if gjk is empty. Since
(8.14) must hold for g = gjk (if non-empty), we can write
r Wr [agj /(1 gj )] + [agk /(1 gk )] [agjk /(1 gjk )]. (8.19)
rGjk
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
which violates (8.14) for g = Gjk . Thus we must have (perhaps after
renumbering)
g1 = {i1 }
g2 = {i1 , i2 }
gR1 = {i1 , i2 , . . . , iR1 }
gR = {i1 , i2 , . . . , iR } = {1, 2, . . . , R}
Proof of Lemma 8.2. This proof is almost identical to the above and need
not be given in full. It suces to note that the right-hand side of (8.15) is
also of the form
br 1 r
rg rg
H1 H1 (8.21)
and
H2 H2 . (8.22)
H1 = H1 = H1 . (8.23)
Similarly, if it can be shown that the set H2 is convex, it would follow that
H2 = H2 = H2 . (8.24)
W = W1 + (1 )W2 ; 0 1. (8.25)
Equations (8.23) and (8.24) are now established; these are important results
which will be referred to as characterisation theorems.
H = {W | S ; S W}.
j 0 (j = 1, 2, . . . , R!), j 0 (j = 0, 1, . . . , R),
R!
R!
0 + j = 1 and +
j Qj = W
j=1 j=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
(using only the rst R1 elements of Qj and W). An initial basis for (8.28)
is obtained by setting j = 0 (j = 1, 2, . . . , R!), 0 = 1, = W.
When an
objective value of 1 is reached (as we know it will be if W is achievable),
the corresponding j S and Qj S dene a mixing strategy S 1 which
achieves W.
Note that there may be (and probably are) many solutions to this
problem. Figure 8.3 shows an example of the set H1 for R = 3 (the priority
ordering at each vertex is indicated in brackets) and a target vector W in
the interior of H1 . This particular target vector can be achieved by mixing
the priority disciplines (123), (321) and (213); or by mixing the disciplines
(213), (132) and (231), etc. The gure suggests that there are eight mixing
strategies which achieve W.
The other case of interest is when the target W is not achievable but
dominates an achievable performance vector (i.e. it satises the inequalities
(8.14) but lies above the conservation law hyperplane). Then, before solving
the linear program (8.28), we have to nd an achievable performance
vector W dominated by W. This is another initial basis problem: W
has to satisfy the inequalities (8.14) plus the conservation law equality,
plus the inequality W W. Again one can introduce articial variables
Fig. 8.3.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
and solve an auxiliary linear program. Since the solution of that program
is a vertex of the set of feasible vectors, the vector W thus obtained is
either one of the vectors Qj (j = 1, 2, . . . , R!), or it satises W r = W r for
some r = 1, 2, . . . , R. This suggests the following alternative algorithm for
nding W.
For r = 1, 2, . . . , R check whether the projection of W along the r-th
coordinate axis into the conservation law hyperplane is achievable; if yes,
then take that projection as W and stop. For j = 1, 2, . . . , R! check whether
Qj is dominated by W; if yes, take Qj as W and stop.
That algorithm may, in some cases, be more ecient than the linear
programming one. For example, if one of the elements of W is obviously
too large, a projection along the corresponding coordinate axis is likely to
yield the result.
The above results apply, with straightforward modications, to M/G/1
systems where pre-emptions are disallowed. The family 2 of scheduling
strategies formed by mixing up to R of the R! head-of-the-line priority
disciplines, is M/G/1-complete. To nd a strategy from 2 which achieves
a pre-specied (and achievable) performance vector W, one solves a linear
program similar to (8.28); the vectors Qj (j = 1, 2, . . . , R!) are replaced
by the performance vectors of the head-of-the-line priority disciplines. If
the target W is not achievable but dominates an achievable performance
vector, one such vector can be obtained either by solving an initial basis
problem or by a vertex searching algorithm.
So, the families 1 and 2 have several attractive features: they are
conceptually simple, easily implementable, parametrised, complete; there
are algorithms for selecting a strategy that achieves (or improves upon) a
given performance vector. We should point out, however, that these mixing
strategies have one important disadvantage: the variances in response times
which they introduce may be unacceptable large, especially in heavily
loaded systems. Suppose, for example, that in an M/M/1 system with two
job classes, the pre-emptive priority disciplines (1, 2) and (2, 1) are mixed in
the proportion = 0.9. Suppose, further, that the system is heavily loaded,
most of the load being contributed by class 2 (say 1 = 0.15, 2 = 0.8).
Then, while most of the class 1 jobs have very short waiting times, a
signicant proportion (approximately 10%) will have to wait much longer;
the over-all mean response time may be as required, but the variance will
be rather large. Managers of computer installations usually avoid such
strategies because the unlucky 10% of the users tend to be more vociferous
than the satised 90%.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
A = lim
0
A,E
E
and therefore
H = lim
0
H,E .
E
B = lim
0
B,E .
E
The proof can now be completed by remarking that, since the surfaces
B,E approach the boundary of H1 , every performance vector W in the
interior of H1 is in the interior of some B,E . As we have seen, this implies
that W H,E and hence W H . Thus H contains all points in H1
except its boundary.
R
j
Wr (t) =1+ nj (t), r = 1, 2, . . . , R. (8.31)
j=1
r
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
Finally, substituting the sums of n and n into (8.31) yields a system
of integrodierential equations:
R
j j
Wr (t) = 1+ e j j t/r
ej u Wj (u)du
j=1
r 0
t
j j (tu)/r
+ e Wr (u)du r = 1, 2, . . . , R. (8.32)
0
R
1 j j
Wr = + (Wj + Wr ) , r = 1, 2, . . . , R. (8.33)
r j=1 j j + r r
R
C(W) = cr Wr , cr 0; r = 1, 2, . . . , R. (8.34)
r=1
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
From the characterisation theorems of section 8.3 we know that the set of
achievable performance vectors is a polytope (H1 in the case of M/M/1
systems, H2 for M/G/1 systems). Moreover, we know exactly what are the
vertices of that polytope (Lemmas 8.1 and 8.2). Since the minimum of a
linear function over a polytope is always reached at one of the vertices, we
immediately obtain the following results.
Lemma 8.4. In an M/M/1 system any cost function of the type (8.34) is
minimised by one of the R! pre-emptive priority disciplines.
j;
Wj < W j+1 ;
Wj+1 > W
r
Wr = W for r < j or r > j + 1.
j + j+1 W
j Wj + j+1 Wj+1 = j W j+1 . (8.35)
Suppose that (cj+1 /j+1 ) > (cj /j ). Then, since the bracketed term which
multiplies (cj+1 /j+1 ) is positive,
> cj (j Wj j W
C(W) C(W) j + j+1 Wj+1 j+1 W
j+1 ) = 0
j
(ii) choose the discipline which gives highest priority to class r1 , second
highest priority to class r2 , . . . , lowest priority to class rR .
R
r
C(W) = Wr ,
r=1
1 M2
WFIFO = + (8.37)
2(1 )
where M2 is the second moment of the service time distribution (see (8.11)).
The restriction on scheduling strategies can be avoided by introducing
articial job classes. Assume rst that the service times can take only
a nite number of values: they are equal to xj with probability pj , j =
1, 2, . . . , J; p1 + p2 + + pJ = 1. The system can be treated as an M/G/1
queue with J job classes, class j having arrival rate j = pj , mean service
time (1/j ) = xj , trac intensity j = pj xj and second moment of service
time distribution M2j = x2j . Scheduling strategies which use information
about exact service times are now admissible because that information is
contained in the class identiers. To minimise the over-all average response
time one has to give class j non-pre-emptive priority over class k if xj < xk .
This is the Shortest-Processing-Time-rst discipline (SPT): when the server
is ready to begin a new service it selects the shortest job of those present
in the system.
The above argument generalises easily to an arbitrary service time
distribution. Thus, in any non-pre-emptive M/G/1 system the average
response time is minimised by the SPT discipline.
We see that the optimal strategy to be followed depends on the amount
of information available. If only the distribution of service times is known
then the best strategy is SEPT which, in the case of one job class, reduces to
serving jobs in order of arrival; the resulting average response time is given
by (8.37). If individual service times are known then the optimal strategy
is SPT; the corresponding average response time is (see expression (1.79),
Chapter 1)
1 M2 dF (x)
WSPT = + )][1 G(x+ )]
(8.38)
2 0 [1 G(x
Fig. 8.4.
Lemma 8.6. The SRPT scheduling strategy achieves the lowest possible
value of the over-all average response time in a G/G/1 queueing system.
time t. The average amount of processor time that the job will actually use
is equal to
y
Qr (t, y) = [1 Fr (t + x)]dx [1 Fr (t)],
0
where Fr (x) is the required service time distribution function for class r.
The probability that the job will complete within the allocated time y is
equal to
These two quantities are used to dene the rank vr (t) of a class r job with
attained service t:
Qr (t, y)
vr (t) = min ,
y cr Sr (t, y)
where cr is the cost coecient associated with class r. The minimum is
taken over the set of permissible allocations y (if jobs can be interrupted
at any point, all y 0 are permissible). The smallest value of y for which
the minimum is reached is called the rank quantum. At each scheduling
decision point, the processor is assigned to the job with the smallest rank
for the duration of the corresponding rank quantum.
When vr (t) is a non-increasing function of the attained service t for
every r = 1, 2, . . . , R, the SR strategy behaves like a pre-emptive priority
discipline based on the ordering (8.36), except that the average service times
are replaced by the average remaining service times. This tends to happen
when the service time distributions have coecients of variation not greater
than 1 (e.g. uniform, Erlang or exponential distributions). If, on the other
hand, vr (t) is an increasing function of t for every r, then SR behaves like a
Processor-Sharing discipline (e.g. for hyperexponential distributions). This
conrms the intuitive idea that when the variations in service times are
small, jobs should not be interrupted too much, and when the variations
are large, it is better to interrupt them often.
The optimality of the SR discipline for the case of no arrivals is proved
by induction on the number of jobs present. We shall omit that proof here;
the interested reader is referred to [13].
A few words should be said about minimising cost functions which are
non-linear in the elements of W. The general problem (min C(W), subject
to the constraints W H1 , or W H2 , depending on whether the system
is M/M/1 or M/G/1) can be tackled by classic mathematical programming
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch08
methods. The fact that the constraints are always linear they are given
by Theorems 8.4 and 8.5 may facilitate the solution. Take, as an example,
a two-class M/M/1 system and consider the problem
min(W12 + W22 )
(nd the point (W1 , W2 ) which is closest to the origin). The constraints are
1 = (K/2)1 ;
W 2 = (K/2)2 ,
W
References
1. Coman, E. G., Jr. and Mitrani, I. (1980). A characterisation of waiting time
performance realisable by single-server queues. Operations Research.
2. Fayolle, G., Iasnogorodski, R. and Mitrani, I. (1978). On the sharing of a
processor among many job classes. Research Report, No. 275, IRIA-Laboria;
also to appear (1980), JACM.
3. Fife, D. W. (1965). Scheduling with random arrivals and linear loss functions.
Man. Sci., 11(3), 429437.
4. Kleinrock, L. (1965). A conservation law for a wide class of queueing
disciplines. Nav. Res. Log. Quart., 12, 181192.
5. Kleinrock, L. (1967). Time-shared systems: A theoretical treatment.
J.A.C.M., 14(2), 242261.
6. Kleinrock, L. (1976). Queueing Systems. Vol. 2. John Wiley, New York.
7. Kleinrock, L. and Coman, E. G., Jr. (1967). Distribution of attained service
in time-shared systems. J. Comp. Sys. Sci., 287298.
8. Mitrani, I. and Hine, J. H. (1977). Complete parameterised families of job
scheduling strategies. Acta Informatica, 8, 6173.
9. ODonovan, T. M. (1974). Distribution of attained and residual service in
general queueing systems. Operations Research, 22, 570575.
10. Schrage, L. E. (1968). A proof of the optimality of the SRPT discipline.
Operations Research, 16(3), 687690.
11. Schrage, L. (1970). An alternative proof of a conservation law for the queue
G/G/1. Operations Research, 18, 185187.
12. Schrage, L. E. and Miller, L. W. (1966). The queue M/G/1 with the shortest
remaining processing time discipline. Operations Research, 14, 670683.
13. Sevcik, K. C. (1974). A proof of the optimality of smallest rank scheduling.
J.A.C.M., 21, 6675.
14. Smith, W. E. (1956). Various optimisers for single-stage production. Nav.
Res. Log. Quart., 3, 5966.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
Chapter 9
Control of Performance
in Multiple-Resource Systems
269
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
Fig. 9.1.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
2b
e(m) = ; b, c > 0 (9.2)
1 + (c/m)2
(Chamberlin, Fuller and Lin [4]). The two functions are shown in Fig. 9.2.
When a new job is admitted into the inner system, i.e. becomes active,
it is allocated a certain number of page frames. This implies, in general, that
the allocation of one or more other jobs is reduced; they begin, therefore, to
operate at a dierent point on their lifetime curves closer to the origin
and start visiting node 1 more often. Similarly, when a job departs from
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
Fig. 9.2.
the inner system, one or more other jobs start visiting node 1 less often.
Thus, not only the service times at node 0 but also the probability that a
job goes to node j after leaving node 0 depend on the number and types of
active jobs.
This behaviour of the inner system means, unfortunately, that the
results of Chapter 3 are not directly applicable to the present model. No
matter what assumptions we make about queueing disciplines, required
service time distributions, etc., a queueing network model of a multi-
programmed system such as the one in Fig. 9.1 will not have a product-
form solution; and without a product-form solution we cannot hope to
have ecient methods for the exact evaluation of performance measures. It
is almost imperative, therefore, to look for approximate solutions.
The parameters of most real-life systems are such that the corre-
sponding models lend themselves easily to decomposition. The interactions
between the inner and the outer systems are weak compared to those within
the inner systems: jobs are admitted into, and depart from, the inner system
at a much lower rate than that at which they circulate inside it. This
allows one to assume that the inner system reaches equilibrium in between
consecutive changes in the degree of multiprogramming. For a given set of
active jobs, the inner system can be treated as a closed queueing network;
that network can be analysed in the steady-state to obtain the rates at
which jobs of various classes obtain CPU service; those rates can be used
to replace the whole inner system by a single server whose rate of service
depends on the system state (via the set of active jobs).
We shall elaborate further on this approach when we apply it to specic
synthesis problems. Our interest will be directed primarily towards the
design of admission and memory allocation policies.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
Fig. 9.3.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
Geometrically, m is such that the ray from the origin passing through the
point (m , e(m )), dominates the entire curve {e(m), m 0}. For example,
the knee of the lifetime function dened by (9.2) is at m = c (see Fig. 9.2).
In order to give the intuitive justication for the knee criterion we
need to introduce a quantity called space-time product, and establish a
relation between it and the . The space-time product, Y , for a job whose
execution time is D (this is total real time spent in the inner system, not
virtual CPU time), is dened as
D
Y = m(t)dt (9.5)
0
where m(t) is the number of page frames that the job holds at time t of its
execution. If m
is the average number of page frames held by the job, then
Y = mD.
(9.6)
MV. On the other hand, the average number of jobs executed during the
period is T (n)V (since an average of T (n) jobs depart per unit time). Hence,
the average space-time product per job is
MV M
Y = = . (9.7)
T (n)V T (n)
1
D= + + (9.8)
e(m)
1 m
Y = m( + 1) + . (9.9)
e(m)
dY /dm = 0, we obtain
1
mopt = c/(1 + 2b( + 1)/ ) 2 .
Hence if 2b( + 1)/ is not large (usually it is less than 1), mopt c = m ;
the optimal allocation is indeed close to the knee of the lifetime function.
An implementation of the knee criterion would involve, for each active
job, a continuing monitoring of its paging activity, an estimation of its
lifetime function and an allocation of memory corresponding to the knee
point. The degree of multiprogramming is thus controlled indirectly via the
memory allocation. Such a control policy can be expected to be expensive,
both in instrumentation and in overheads. On the other hand, as Denning
and Kahns experiments suggest (see [6]), the knee criterion is robust and
yields near-optimal degrees of multiprogramming over a wide range of
loading conditions. One cannot, of course, apply the knee criterion if jobs
behave according to a lifetime function that has no knee, such as the one
dened by (9.1). A nite value for mopt may still exist, as we shall see
shortly, and it may be possible to obtain an estimate for it. The memory
allocation controller can then use that estimate.
where ej is the j-th most recent interval, regardless of which job gener-
ated it.
The L = S control policy attempts to balance the system lifetime
and the paging drum average service time S: it maintains the degree of
multiprogramming n at such a level that
where c is a constant not much greater than 1. The intuition behind this
criterion derives from the bounds that device service times place on .
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
Let 1/bi and Ui be, respectively, the mean service time and the
utilisation of node i in the inner system (i = 0, 1, . . . , K; 1/b1 = S).
The parameters b1 , b2 , . . . , bk are independent of n but b0 depends on it via
the memory allocation and lifetime functions. The utilisations Ui depend, of
course, on n. Denote further by 1/ai the mean CPU interval (not necessarily
continuous) between consecutive requests for device i; i = 1, 2, . . . , K. These
are averages over all active jobs. 1/a1 is the system lifetime L(n); all other
ai s are program characteristics and are independent of n.
While the CPU is busy, jobs depart from it in the direction of node i
at rate ai (i = 1, 2, . . . , K); the probability that the CPU is busy is
U0 ; therefore, the rate at which jobs arrive into device i is equal to
U0 ai . Similarly, the rate at which jobs leave device i is equal to Ui bi
(i = 1, 2, . . . , K). Since these two rates are equal in the steady-state, we have
U 0 ai = U i b i i = 1, 2, . . . , K. (9.12)
From the above, the is bound by two functions. One of them is constant
and the other is decreasing with n (the more active jobs, the less memory
per job and hence the smaller CPU intervals between page faults). These
bounds are illustrated in Fig. 9.4. If the two bounds intersect, they do so at
point n
which satises L( n) = IS. Intuitively, the should start decreasing
for n > n , i.e. n
is slightly larger but close to the optimal value of n.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
Fig. 9.4.
(k 1)
amkopt = .
+1
If the system is not I/O-bound, the constant (the average I/O delay per
unit of CPU time) is small; the constant k of the Belady lifetime function
is usually less than 2. Assuming that (the average delay per page fault) is
not much greater than the drum service time S, we conclude that e(mopt )
is near S: we arrive again at the L = S criterion.
A controller based on the L = S rule is simpler and easier to implement
than one based on the knee criterion. All that is required is an estimate
of the current system lifetime L(n), obtained as in (9.10). However, this
policy appears to be less robust and in some cases (especially in I/O-bound
systems) leads to a degree of multiprogramming which is signicantly lower
than the optimal (see [6]).
The supporting argument for the 50% criterion is equally simple. It proceeds
as follows.
One of the manifestations of thrashing is a high rate of page faults,
hence a high rate of requests for the paging drum. A high rate of requests
implies a long queue; thus, a queue at the paging drum is symptomatic of a
thrashing system and, to prevent thrashing, queues should not be allowed
to develop. An average of one request at the paging drum is a good target
to aim for. If we treat that device as an independent M/M/1 queue we can
write for the average number n of requests there (see Chapter 1)
= U1 /(1 U1 ),
n
where U1 is the trac intensity, or the probability that the drum is busy. If
we wish that average to be n = 1 we should keep the utilisation at U1 = 0.5.
Tenuous though the above argument may appear, the 50% rule seems
to perform reasonably well, especially in systems where the average drum
service time S is lower than the knee lifetime e(m ) (see [6]). In other
cases this criterion tends to be less reliable than the other two and to
underestimate the optimal degree of multiprogramming. On the other hand,
it is the simplest of the three and the most straightforward to implement.
To summarise, we have described here three heuristic control rules
for maintaining the degree of multiprogramming at, or near, its optimal
level. We say that they are heuristic because there are no mathematical
proofs establishing their validity, only intuitive arguments. There is,
however, a certain amount of empirical and numerical evidence [1, 12, 13]
which suggests that these rules can be applied successfully in practical
systems.
It should be emphasised that any dynamic control is necessarily
involved with transient phenomena whereas the theoretical support for the
proposed control procedures is based on steady-state analysis. The degree
of multiprogramming is assumed to remain constant long enough for the
inner system to reach steady-state. If the loading conditions change rapidly
in relation to the control actions (or, alternatively, if the control procedure
reacts slowly to changes in the load conditions) then this assumption is
violated and the control, if it works at all, will be unstable. There are no
general assertions that can be made in this connection; much depends on
implementation and instrumentation, as well as on the control algorithm.
For example, some experiments (Leroudier and Potier [11]) indicate that
the 50% rule responds rapidly to load uctuations.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
Fig. 9.5.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
jobs should remain impeded longer; when the improves, impeded jobs can
be reintroduced into the active set at a higher rate. Under the page fault
rate control policy, the average time 1/ that jobs remain in the impeded
set is proportional to the average interval 1/ between departures from the
active set. We set
= h (9.17)
where h is a small constant (in the numerical evaluations of the policy its
value was chosen as h = 0.01).
Thus, an implementation of the policy would involve two types of
monitoring: (i) counting page faults for each active job in order to decide
when to remove jobs to the impeded set, and (ii) estimating the total rate
of departures from the inner system in order to regulate the average times
that jobs spend in the impeded set. This compromise between program-
driven and load-driven control allows the policy to prevent thrashing by
discriminating against the jobs which contribute to it most the jobs
with most page faults.
An exact analysis of the system performance under RCP is not feasible
for the reasons mentioned in the last section: the behaviour of active jobs
depends on their number. However, an approximate evaluation can be
obtained by applying decomposition. First, we consider the inner system
as a closed network with a xed number nr of class r jobs circulating
inside (r = 1, 2, . . . , R). An analysis of the closed network will enable
us to replace the whole inner system by a single aggregate server which
gives simultaneous service to all active jobs, at rates depending on the
state n = (n1 , n2 , . . . , nR ). To do this we need the steady-state probability
r (n1 , n2 , . . . , nR ) that, in the closed network, a class r job is in service
at the CPU: that probability determines the for class r jobs in state n.
Also necessary is the steady-state probability rJ (n1 , n2 , . . . , nR ) that a
class r job which has already had J page faults is in service at the CPU: it
determines the rate at which jobs leave the inner system to join the impeded
set.
An important distinguishing feature of the jobs of a given job class is
their paging behaviour. In the model, a dierent lifetime function er (m)
is associated with each job class (r = 1, 2, . . . , R). The counting of page
faults in the closed network is modelled by splitting class r into J + 1
articial job classes (r, 0), (r, 1), . . . , (r, J). At each visit to the paging
drum, a job of class (r, j) becomes a job of class (r, j +1) if j = 0, 1, . . . , J 1
and it becomes a job of class (r, 0) if j = J.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
Exact expressions for r (n) and rJ (n) can be obtained, under a suitable
set of assumptions, by applying the BCMP theorem of Chapter 3 (the
formulae for a special case can be found in [8]). However, the computational
eort associated with the solution of a multiclass network (calculation of
normalisation constant, aggregation of states to nd marginal distributions,
etc.) is considerable. What is more, that eort grows rather rapidly with
the size of the model, in particular with the number of articial job classes
resulting from a large value of J. On the other hand, an exact solution
is rarely necessary, especially in view of the fact that the whole model is
approximate (because of the decomposition and because parameters have to
be estimated). Good approximations for the probabilities r (n) and rJ (n)
can be obtained rather easily as follows.
Solve a single-class closed network with n = n1 + n2 + + nR jobs
circulating inside. As a lifetime function, use the linear combination
R
e= (nr /n)er .
r=1
If other parameters vary across job classes, they are also averaged in a
similar fashion. That single-class solution yields the over-all CPU utilisation
U0 (n).
Now return for a moment to the full model (including the outer system
and the impeded set) and let 1/r be the average total CPU time required
by a class r job; assume that the distribution of that time is exponential.
Then, while a class r job is being served by the CPU, it leaves for the outer
system at rate r . Similarly, while a class r job is in service at the CPU it
leaves for the impeded set at rate 1/[(J + 1)er ] (it is ejected at the J + 1st
page fault). Therefore, in the virtual CPU time of a class r job the average
interval between consecutive departures from the inner system is equal to
1
1 (J + 1)er
Qr = r + = ; r = 1, 2, . . . , R. (9.18)
(J + 1)er 1 + (J + 1)er r
Intuitively, the proportion of all CPU busy time which is devoted to
class r jobs can be approximated by
R
qr = nr Qr ns Qs ; r = 1, 2, . . . , R. (9.19)
s=1
Furthermore, the class r jobs which have had J page faults (these are
the class (r, J) jobs) occupy approximately a fraction 1/(J + 1) of that
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
r (n1 , n2 , . . . , nR ) = qr U0 (n1 , n2 , . . . , nR )
rJ (n1 , n2 , . . . , nR ) = [qr /(J + 1)]U0 (n1 , n2 , . . . , nR ); r = 1, 2, . . . , R.
(9.20)
The total rate at which class r jobs depart from the aggregate server is
r (n) = r (n) + r (n). The rate at which each impeded job returns to the
aggregate server is given by (9.17), with = 1 + 2 + + R .
The model has thus been reduced to a queueing network with three
nodes: the outer system, the aggregate server and the impeded set. The
state of the network (under exponential assumptions) is a 2R-dimensional
Markov process (n; k) = (n1 , n2 , . . . , nR ; k1 , k2 , . . . , kR ), where nr and kr
are the numbers of class r jobs at the aggregate server and in the impeded
set, respectively (r = 1, 2, . . . , R). While it is easy to write the steady-state
balance equations for that process, solving them is by no means easy: there
is no product-form solution because the behaviour of the aggregate server
depends on the vector n and not just on the total number of jobs there.
However, in some cases another level of decomposition may simplify matters
considerably.
Consider the case of two job classes and suppose that the outer system
consists of two nite sources: N1 terminals of class 1 and N2 terminals
of class 2. There are thus N1 and N2 class 1 and class 2 jobs circulating
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
endlessly among the three nodes; the system state is (n1 , n2 ; k1 , k2 ), where
n1 + k1 N1 , n2 + k2 N2 . The number of states, and hence the number
N1 + 2 N2 + 2
of balance equations, is 2 2 .
If the two job classes have very dierent lifetime function character-
istics, we can attempt another decomposition. Suppose that the ejection
threshold J can be chosen in such a way that jobs of one class, say class 2,
are ejected from the inner system much more often than the others (the
choice of J will be examined later). Then, if the aggregate server and the
impeded set are considered in isolation, with I1 and I2 class 1 and class
2 jobs circulating among them, it can be assumed that all I1 jobs are at
the aggregate server. The probabilities P (I1 , n2 | I1 , I2 ) that there are n2
class 2 jobs at the aggregate server, given I1 and I2 (n2 = 0, 1, . . . , I2 ) can
be obtained by solving a simple I2 + 1 state Markov process.
As the next step, the aggregate server and the impeded set are replaced
by a single server which, when in state I1 , I2 , sends class r jobs to the
terminals at rate
I2
r (I1 , I2 ) = r (I1 , n2 )P (I1 , n2 | I1 , I2 ), r = 1, 2,
n2 =0
Fig. 9.6.
parameter k of the lifetime function: k = 1.8 for the good class 1 jobs
and k = 1.5 for the bad class 2 jobs. Main memory was divided equally
among the active jobs and an ejection threshold J = 30 was used.
On the basis of some numerical comparisons (see [7] and [8]) it appears
that if the load characteristics do not vary rapidly with time, the control
exercised by the page fault rate control policy is close to optimal. Those
evaluations, however, ignored the overheads associated with the policy and,
in particular, the times taken to move jobs between the active and impeded
sets. The ability of the policy to react quickly to variations in program
behaviour is also open to investigation.
The choice of J: The function of the control parameter J is twofold:
rst, to prevent thrashing eectively and second, to ensure that the good
jobs those that do not have page faults often are not ejected from
the active set often. To full the rst objective J should not be too large
(otherwise there would be no control), and to full the second it should not
be too small. This trade-o can be assessed by evaluating the probability
pr that a class r job is ejected from the active set before it is completed.
Under exponential assumptions about total execution and inter-page fault
times we can write (see Chapter 1)
J J
1/er 1
pr = = , r = 1, 2, . . . , R.
r + (1/er ) er r + 1
In the case of two job classes suppose that class 1 is much better than
class 2 (i.e. e1 e2 ) and that it is possible to choose J so that Je1 1 1
and Je2 2 1. Then we would have p1 0, p2 1 and both objectives
would be satised. If the dierence between the two classes is not so extreme
and it is not possible to achieve p1 0 without at the same time having
p2 0, then the decision is much less clear-cut. One could proceed by
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
The memory allocated to class r is divided equally among the active class r
jobs. This memory allocation strategy must be accompanied by a job
admission strategy, in order to avoid thrashing.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
One could then decide to admit a class r job into the active set if there are
cr free pages in the allocation for class r (the knee criterion, section 9.3).
In this last case, the number of class r jobs in the active set would be
Here, as before, 1/r is the average CPU time required by a class r job;
the distribution of that time is assumed to be exponential.
The global system state is now described by the vector N =
(N1 , N2 , . . . , NR ), where Nr is the number of class r jobs submitted for
execution. Using the appropriate mapping N n (equation (9.25) is an
example) in conjunction with (9.26), one can write balance equations for the
steady-state distribution of N. These equations have to be solved numeri-
cally (perhaps using an approximation technique: see [10]), since there are
no closed-form solutions for R-dimensional Markov processes. From the dis-
tribution of N one can nd the mean number of class r jobs submitted and
hence, by Littles theorem, the average response time Wr (r = 1, 2, . . . , R).
Some results for the case of two-job classes are illustrated in Fig. 9.7.
The outer system in that example consisted of two independent Poisson
streams of jobs, with rates 1 and 2 for class 1 and class 2, respectively.
The inner system comprised a CPU, a paging drum and a ling disk;
Fig. 9.7.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
Fig. 9.8.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
The total CPU utilisation U0 is obtained by summing (9.28) over all job
classes:
R
Nr s0r
U0 = . (9.29)
r=1
Wr + r
When U0 tends to zero, that curve moves away from the origin and attens
out; conversely, when U0 increases the curve moves towards the origin and
becomes more convex.
Thus, in a well-tuned system where the CPU utilisation is close to the
maximum attainable, the performance vectors of all scheduling strategies
which maintain that utilisation lie on a hyperbola. Of course, not all points
on the hyperbola are achievable. As in the single-server case, the achievable
performance vectors must satisfy inequalities of the type
Wr Wrmin , r = 1, 2,
where Wrmin is the average response time for class r jobs in a system where
the other job class does not exist and where the CPU utilisation is the
maximum attainable. The set of achievable performance vectors is therefore
contained in a region such as the shaded area in Fig. 9.9. The situation is
similar when the number of job classes is R > 2.
We should point out that this characterisation of achievable perfor-
mance has a limited, if any, practical value. The maximum attainable CPU
utilisation is not usually known and the position of the bounding hyperbola
(or R-dimensional surface) may be very sensitive to its estimate. Moreover,
even if the constraints are calculated accurately there is no guarantee that
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
Fig. 9.9.
all points which satisfy them are, in fact, achievable. What we have obtained
is an idea of the likely shape of the set of achievable performance vectors.
That idea is consistent with the results of the last section (the curve in
Fig. 9.7b closely resembles a hyperbola). It also conrms the observation
made before, namely that the trade-os between response times for dierent
job classes are likely to be non-linear.
References
1. Adams, M. C. and Millard, G. E. (1975). Performance Measurements on
the Edinburgh Multi-Access System (EMAS). Proc. ICS 75, Antibes.
2. Belady, L. A. and Kuehner, C. J. (1969). Dynamic space sharing in computer
systems. Comm. A.C.M., 12, 282288.
3. Buzen, J. P. (1976). Fundamental operational laws of computer system
performance. Acta Informatica, 7, 167182.
4. Chamberlin, D. D., Fuller, S. H. and Lin, L. Y. (1973). A Page Allocation
Strategy for Multiprogramming Systems with Virtual Memory. Proc. 4th
Symp. on Operations Systems Principles, pp. 6672.
5. Denning, P. J. and Kahn, K. C. (1975). A Study of Program Locality and
Lifetime Functions. Proc. 5th Symp. on Operations Systems Principles,
pp. 207216.
6. Denning, P. J., Kahn, K. C., Leroudier, J., Potier, D. and Suri, R. (1976).
Optimal multiprogramming. Acta Informatica, 7, 197216.
7. Gelenbe, E. and Kurinckx, A. (1978). Random injection control of multi-
programming in virtual memory. IEEE Trans. on Software Engng., 4, 217.
8. Gelenbe, E., Kurinckx, A. and Mitrani, I. (1978). The Rate Control Policy
for Virtual Memory Management. Proc. 2nd Int. Symp. on Operations
Systems, IRIA, Rocquencourt.
9. Hine, J. H. (1978). Scheduling for Pre-specified Performance in Multipro-
grammed Computer Systems. Research. Report., University of Wellington.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
10. Hine, J. H., Mitrani, I. and Tsur, S. (1979). Control of response times in
multi-class systems by memory allocation. C.A.C.M., 22(7), 415423.
11. Leroudier, J. and Potier, D. (1976). Principles of Optimality for Multipro-
gramming. Proc. Int. Symp. on Computer Performance Modelling, Measur-
ing and Evaluation, pp. 211218. Cambridge, Massachusetts.
12. Rodriguez-Rossel, J. and Dupuy, J. P. (1972). The Evaluation of a Time
Sharing Page Demand System. Proc. AFIPS, SJCC 40, pp. 759765.
13. Sekino, A. (1972). Performance Evaluation of a Multiprogrammed
Time Shared Computer System. MIT Project MAC, Research Report
MAC-TR-103.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch09
Chapter 10
A Queue with Server of Walking Type
10.1. Introduction
Queues with autonomous service (QAS) represent service systems in which
the server becomes unavailable for a random time after each service epoch.
Such systems have been used to model secondary memory devices in
computer systems (e.g. paging disks or drums) as was done in Chapter 2.
The queue with server of walking type studied by Skinner [1] is a special
instance of our model. This model has also been considered by Borovkov [4].
Assuming general independent interarrival times we obtain an opera-
tional formula relating the waiting time in stationary state of a QAS to
the waiting time of the GI/G/1 queue. This result dispenses the need for
analysis of the QAS in special cases and generalizes the result of Skinner [1],
or that of Coman [2] for a paging drum. Sucient conditions for stability
or instability of the system are also obtained.
297
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch10
That is, service will resume for the (k + 1)-th customer which arrives at
time ak+1 k+1 1 Ai at time sk + Tk + 1l(ak+1 ) Sik , where
l
l(ak+1 ) = inf l : sk + T k + Sik ak+1 .
1
We assume that the {Snk }n,k1 are i.i.d. and independent of the interarrival
times and of the sequence {Sn }n1 . In the sequel, we shall drop the index
k associated with Snk in order to simplify the notation, though it will be
understood that the variables associated with the end of dierent busy
periods are distinct.
The model we consider arises in many applications. In computer
systems [2, 3, 5] it serves as a model of a paging drum (in this case S
and S are constant and equal). In data communication systems it can
serve to represent a data transmission facility where transmission begins
at predetermined instants of time.
Using the terminology of Skinner [1] who analyzed the model assuming
Poisson arrivals, we shall call it a queue with server of walking type: after
each service the server takes a walk. Borovkov [4] studies a related model
which he calls a queue with autonomous service.
The purpose of this paper is to obtain a general formula relating the
waiting time Wn of the n-th customer in our model to the waiting time,
of the n-th customer Vn in an equivalent GI/G/1 queue, n 1. This
equivalent GI/G/1 queue has the same arrival process, but the service times
are S1 , S2 , . . . , Sn , . . . and Vn+1 = [Vn + Sn An+1 ]+ . This result allows
us to dispense with a special analysis of our queueing model in stationnary
state since we can obtain the result directly from the known analysis of the
corresponding GI/G/1 queue.
The formula (Theorem 4) is derived in section 2 together with sucient
conditions for ergodicity. Section 3 contains an application to the paging
drum model.
distribution of the Sn and the An . His main result is that the queue
length distribution (where the queue does not include the customers in
service) of the above system is identical to the queue length distribution of
a conventional queue (with batch arrivals and batch service) if the service
times are exponentially distributed. The model considered by Skinner [1] is
a special case of the one we study since he assumes that the arrival process is
Poisson; otherwise it is identical to ours. He obtains the generating function
for the queue length distribution in stationary state.
Lemma 1.
Proof. Notice from (10.1) that (x) x for all x with probability 1: if
x 0 the statement is obvious; since (x) 0 with probability 1 it follows
that (x) x if x > 0. Therefore, by Lemma 1 we have
Wn+1 Wn + n , n1
n
Therefore Wn+1 1 n , n 1. If En > 0, then the sum on the RHS
converges with probability 1 to + as n .
Henceforth we shall assume that En < 0 for all n 1.
We shall now study the characteristic function EeitWn+1 for the waiting
time process. Using (10.1) we have, for any real t
Pk+1
EeitWn+1 = Eeit(Wn +n ) I[Wn + n 0] + Eeit(Wn +n + 1 Si )
.I
k=0
k+1 k
Wn + n + Si 0, 0 > Wn + n + Si (10.5)
1 1
Let f (t) = EeitS .
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch10
Then
Pk
EeitWn+1 = Eeit(Wn +n ) + [f (t) 1] Eeit(Wn +n + 1 Si ) .I
k=0
k
Wn + n + Si < 0 (10.6)
1
(a) the random variable is not arithmetic; that is g(t) = Eeit has a single
real value t(t = 0) for which g(t) = 1,
(b) E < 0, and E S < .
Then:
l(x)
= lim Si x,
p x
1
and is independent of V .
Proof. Dene
k
P
it(Wn +n + k
1 Si )
n (t) = Ee .I Wn + n + Si < 0
k=0 1
k
0
itx
= e d P Wn + n + Si < x (10.7)
k=0 1
n (t) = EeitWn
Our proof will be complete if we can prove the existence and uniqueness
of the characteristic function (t) EeitW of a positive random variable
W , which is the solution of the stationary equation
obtained from (10.8), such that (i) and (ii) are satised.
Uniqueness. We shall rst show that if the solution (t) to (10.9) exists,
then it is unique. If (t) exists, it must be continuous for real t and (t)
must exist. Using (10.7):
0
(t) = eity dG(y)
where
k
G(y) P W ++ Si < y
k=0 1
P x++ Si < 0 < P + Si < 0 = EH()
k=0 1 k=0 1
P W ++ Si < y P +
Si < 0
1 1
Consider the LHS of (10.10). (t), f (t) and EeitV are characteristic
functions of positive random variables; they are therefore analytic in the
upper half-plane (Im(t) > 0) and continuous on the real line and bounded.
Consider the RHS of (10.10). (t) is the characteristic function of a negative
random variable and so is EeitX ; therefore the RHS of (10.10) is analytic
on the lower half-plane (Im(t) < 0) and continuous on the real line and
bounded. Therefore by the Liouvilles Theorem the expression (10.10) is a
constant, call it C. Let us write:
C(f (t) 1) itV
(t) = Ee
itP (V = 0)
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch10
Taking
CE S
1 = (0) =
P (V = 0)
we have C = P (V = 0)/E S and
f (t) 1 itV
(t) = Ee (10.11)
itE S
Therefore if (t) exists, then it is unique since it is given by (10.11).
In fact, we have also shown that if it exists, it satises (ii) since (10.11) is
simply the Fourier transform of the statement in (ii).
Existence. We must now prove the existence of the solution (t) given by
(10.11), of the equation (10.9).
Using (10.7), we shall show that (t) of (10.11) is a solution to (10.9).
We write, from (10.7):
0
k
(t) = d P W ++
Sn < x eitx (10.12)
k=0 1
k
2 [0, x[= P W+ Sn < x .
k=0 n=1
(t) EeitV
=
1 f (t)
itE(S)
We deduce = + 2 , where
2 is restricted to R+ .
is the restriction of to R and therefore has the Fourier transform
(t) which is
1 E(eitX )
(t) = p(V = 0)
itE S
Hence replacing (t) above and (10.11) in (10.9) we see that the equality
(10.9) is satised completing the existence proof.
We have established the existence and uniqueness of the stationary
solution (t) of equation (10.8). We now have to prove that
lim Wn = W
n p
i.e. that this stationary solution is the limit in the above equation. For this
we shall call upon general results on the ergodicity of Markov chains as
presented by Revuz [6]. In particular:
1 We rst show that Wn is irreducible.
2 We use the Theorem (Revuz [6], Theorem 2.7, Chapter 3) that states
that if a chain is irreducible and if a nite invariant measure exists, then
it is recurrent in the sense of Harris (i.e. a Harris chain). Thus we show
that Wn is a Harris chain.
3 Finally we use Oreys theorem (Revuz [6], Theorem 2.8, Chapter 6)
which states that if a nite invariant measure m exists for an aperiodic
Harris chain Wn , then Wn W; if the measure m is a probability
p
measure then it is the measure of W .
1 To show irreducibility, consider the measure m whose Fourier transform
is (t). By (10.11) we can write
m=vs
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch10
for a subset A of the non-negative real line. We shall show that, for each
initial state x [0, [, there exists a positive integer n such that
P (Wn A | W0 = x) > 0.
But since P (Vm = Wm A) > 0 for each nite m (the case where the
queue with automous server does not empty up to, and including, the m-th
customer), then
2 Theorem 2.7, Chapter 3 of Revuz [6] states that Wn is a Harris chain
if it is v-irreducible and if there exists nite invariant measure m such
that v(A) > 0 m(A) > 0 for all A (m v, in Revuzs notation).
This has already been proved. Therefore Wn is indeed recurrent in the
sense of Harris.
3 We now have to show, in order to use Oreys theorem, that Wn is
aperiodic. We call again upon the classical result that Vn is ergodic
if E < 0 and is not arithmetic (both of which we have assumed).
Therefore Vn is aperiodic, and so is Wn since for each nite m
P (Vm = Wm A) > 0
lim Wn = W
n p
W =V + Y
p
References
1. Skinner, C. E. (1967). A priority queueing model with server walking type,
Operations Research, 15, 278285.
2. Coman, E. G. (1969). Analysis of a drum input-output queue under scheduled
operation, J.A.C.M., 16(1), 7390.
3. Gelenbe, E., Lenfant, J. and Potjer, D. (1975). Response time of a xed-head
disk to transfer of variable length, S.I.A.M.J. on Computing, 4(4), 461473.
January 11, 2010 12:17 spi-b749 9in x 6in b749-ch10
Index
309
January 11, 2010 12:18 spi-b749 9in x 6in b749-Index
Index 311
little theorem, 87, 88, 108, 112, 236, node, 7377, 7990, 9395, 98114,
257, 264, 285, 289 193, 195, 196, 199, 271274, 278,
local balance, 85, 92, 93, 97, 102, 105, 284, 285, 292
226 closed, 74, 81, 111
locality of reference, 269 generating function, 107
lost arrival, 106 open, 74, 81, 83, 101
lumpability of stochastic matrices, recurrent, 8183, 87
218 transient, 81, 82
utilisation, 87, 106, 108110,
main memory, 59, 269272, 274, 281, 112, 278, 292
286288, 290, 291 normal distribution, 166, 187
M/D/1 system, 190 normalising constant, 106, 107, 110,
M/G/1 system, 50, 58, 242, 243, 248, 111, 113, 114, 200
249, 252, 260262, 265 normalising equation, 85, 91, 93, 94,
M/G/c system, 80 96, 101
M/M/1 system, 92, 225, 240, 242, Nortons theorem, 224, 226, 229
243, 246, 248250, 252, 253, 260, numerical discretisation, see
261, 265267 Discretisation
M/M/c system, 84
marginal distribution, 110, 212, 213, optimal scheduling strategies, 259
215, 225, 283, 285 output theorem, 78, 79, 90
Markov chain, 43, 45, 49, 50, 82, 86,
page fault, 272, 274, 276, 278,
87, 105, 186, 201, 206, 215, 216
280283, 286
embedded, 43, 45
page fault rate control policy, see rate
Markov process, 85, 96, 97, 100, 284,
control policy
285, 289
page frames, 271, 272, 274276
irreducible, 8587
paging device, 271, 272
transient, 82
paging drum, 44, 58, 6365, 74, 113,
Markov renewal process, 49, 50
275, 277, 279282, 289, 292
Markov renewal theory, 43, 49 path trac, 196
mass service systems, 1 performance evaluation, 1, 58, 201,
measures of system performance, 106 229, 231
memory allocation, see also Selective performance measures, 43, 44, 75, 76,
memory allocation 106, 107, 113, 198, 231, 273
memoryless property, 235 performance objective, 231, 232, 248,
mixing scheduling strategy, see 270
Scheduling strategy performance vector, 232, 236, 237,
MM property, 105 239, 241, 242, 247256, 258260,
multiplexed communication channels, 267, 270, 271, 287, 288, 290294
43, 44 achievable, 231, 232, 236, 237,
multiprogrammed computer system, 239, 242, 247253, 258, 260,
74, 267 270, 288, 290, 291, 293, 294
Poisson
negative customers, 117, 118 arrival stream, 57, 63, 78, 79, 238
network generating function, 109 non-homogeneous, 99
network, see queueing network Poisson process, 59, 90, 99
January 11, 2010 12:18 spi-b749 9in x 6in b749-Index
Index 313