IRIA
laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
78150 - Le Chesnay
Tél. 954.90.20

laboratoire de recherche
en informatique
et automatique

<figure><img src="image_1.png" /></figure>

Probabilistic models
of computer systems
part 1 (exact results)

E. Gelenbe, R.R. Muntz

Rapport de Recherche n°146
Décembre 1975
# PROBABILISTIC MODELS OF COMPUTER SYSTEMS - PART I
(EXACT RESULTS)

E. Gelenbe, R.R. Muntz

<figure><img src="image_2.png" /></figure>

## RÉSUMÉ :

Ce papier donne une vue d'ensemble des résultats importants de la théorie des réseaux de files d'attente. Le choix des sujets présentés a été orienté par les possibilités d'application à l'analyse, par les mathématiques ou la simulation, de systèmes informatiques complexes. Les principaux sujets abordés sont les suivants : les théorèmes de Jackson, et de Gordon et Newell ; les théorèmes de taux d'utilisation de Chang, Lavenberg et Traiger, les résultats de Baskett, Chandy, Muntz, Palacios. Nous rassemblons également les travaux relatifs aux processus de Poisson dans un réseau de file d'attente, pour en faire un ensemble cohérent. Deux articles (part II, III) suivront celui-ci, qui présenteront de nouveaux résultats à l'aide d'approximations par les processus de diffusion, pour étudier les réseaux de files d'attente. Des applications pour l'évaluation des performances de systèmes informatiques y seront données.

## ABSTRACT :

This paper is a survey of certain known results concerning networks of queues. The choice of the topics presented has been made with special emphasis on mathematical results which can be applied to the analysis and synthesis of mathematical or simulation models of complex computer systems in which an ensemble of resources is shared among a set of user programs. The subjects covered include the Jackson, and Gordon and Newell theorems ; the work-rate theorems of Chang, Lavenberg and Traiger ; the Baskett, Chandy, Muntz, Palacios results. We also discuss in a unified manner problems related to Poisson processes in queueing networks. Companion papers (Parts II, III) will present new results concerning diffusion approximations to queueing networks and some of their applications to computer system performance evaluation.
The theory of networks of queues is developing rapidly under the impact of problems which have been raised by the mathematical modelling of multiprogrammed computer systems and related topics. The thrust of recent work has been towards establishing explicit expressions for the equilibrium probabilities of queueing networks of increasing complexity $[1,2,3,4]$; these results serve largely as a guideline to the search for approximate solutions of queueing models for which exact solutions have not yet been found and which may well be unobtainable. Interesting and useful results have been obtained by approximating queueing systems with continuous processes $[5,6,7,8]$ which appear to yield very accurate predictions for wide ranges of parameter values $[9,10]$. Iterative techniques based on heuristic considerations $[11]$ have also been considered and appear quite promising. Approximate solutions to some models in which intricate structural dependencies exist have been given in $[12,13]$. The analysis of systems in which blocking phenomena occur $[16]$ has also been examined.

The purpose of this paper is to present in a unified manner some of the fundamental results which are directly applicable to the performance analysis of multiprogramming systems.

We shall begin with Jackson's theorem [1] concerning a queueing network (QN) with exponentially distributed service time, Poisson arrivals, and routing in the network described by a first-order Markov chain. Both the "closed" and "open" version of the result will be given and the relationship to a representation of computer system behaviour will be presented. The shortcomings of Jackson's theorem for the accurate representation of such systems will motivate the remaining discussion.

In section III the model considered by Jackson is generalized to include arbitrary service time distributions. In this case little can be said concerning equilibrium probability distributions. There is, however, an interesting theorem due to Chang and Lavenberg [17] yielding a relationship between the work-rates of the different processors for a closed QN.

Another work-rate theorem [18] is given stating that work-rates are unchanged if in Jackson's QN some or all of the servers at which queues cannot build up are replaced by servers with service times of general distribution but identical mean.
Section IV is devoted to the presentation of the BCMP (Baskett, Chandy, Muntz, Palacios) [2] theorem.

The result in [2] is a substantial contribution to the area. Jackson's model is extended to one in which several classes of customers share the QN. Let a service station operate under a BCMP discipline if the latter is either a round-robin processor sharing (RR) or a pre-emptive-resume priority last-come-first-served discipline (LCFS). The product form equilibrium solution to the network equations is given [2] for a QN in which almost general distribution functions are allowed at stations containing an infinite number of servers or those utilizing a BCMP discipline.

The fifth section is devoted to an examination of the relationship between Poisson processes and queueing networks. Theorems are presented which show that under some conditions a point process formed by the superposition of several random processes or the sampling of a single random process can approach a Poisson process.

In section 6 we summarize the methods which have been developed for computing utilizations, mean queue lengths, etc., from the equilibrium state probabilities. The large state space of queueing network models requires that these computations be as efficient as possible.

Solution methods for first-come-first-served discipline and general service time distributions are not treated in this paper. Such QN can be approached using diffusion approximations [6,7,8,9]; companion papers will present some new results and applications in this direction.

A survey article by R. Muntz [30] discusses some of the applications which have been made of the methods presented here to the performance evaluation and prediction of existing systems.

## II. JACKSON'S THEOREM

Consider the QN of Figure 1. A set of service stations are interconnected to each other via an arbitrary topology. Customers enter the network
via a source and are directed to a service station. These stations allocate
service to the customers awaiting service in the respective queues in FCFS
order. The analogy with a complex computer system is immediate: the customers
are tasks or processes (in the sense attributed to this word in computer
science) and the service stations are processors (central processing units,
input-output channels, etc...), although a large number of application areas
outside computer systems may be imagined. Another analogy with a computer
system may also be considered: the service stations could be processes
utilizing a set of resources (the customers in Jackson's model); these
resources are created by an external agent (the arrival process) or are de-
stroyed (the departure of customers), or again a resource is obtained by
process j (service station j) after having been utilized by process i.

Let us now turn to a precise definition of the model which consists of
n service stations. Let $k = (k_1, ..., k_n)$ be a vector of non-negative integers,
$k_i$ being the number of customers present at service station i, and denote by
$K = \sum_i^n k_i$ the total number of customers in the network. The customer service
times at station i are assumed to be independent and exponentially distributed
with mean $1/\mu_i$. The rate at which service is given to customers may vary
with the number of customers $k_i$ at station i. Thus the departure rate may
be a function of $k_i$ which will be denoted by $\mu_i(k_i)$.

Let $\gamma_m$ be the m.th service station visited by a customer which has entered the system ; then $\{\gamma_m ; m = 1, 2, ...\}$ is a Markov chain whose transition matrix $P = (p_{ij})$ is defined as $p_{ij} = \text{Pr}[\gamma_{m+1} = j | \gamma_m = i]$ for any $m = 1, 2, ..., $ where i and j correspond to service stations. In the "open" model $p_{oj}$ will be the probability that an external arrival first enters station j if $1 \le j \le n$. We take $p_{o, n+1} = 0$. For $1 \le i \le n$, $p_{i, n+1}$ is the probability that a customer leaves the system after completing service at station i. The arrival process of customers is assumed to be Poisson (interarrival times are exponentially distributed) with rate $\lambda(K)$ i.e. a function of the total number of customers in the system. With these assumptions it is immediately seen that the stochastic process representing the number
of customers in each queue is a continuous time Markov chain [20]. Let $a(k, i)$ and $b(k, i)$ be vectors identical to $k$ except for the element $k_i$ which is replaced by $k_{i+1}$, $k_{i-1}$ respectively. Let also $C(k, i, j)$ be identical to $k$ except for $k_i$ an $k_j$ which are replaced by $k_{i+1}$, $k_{j-1}$ respectively. Define the probability $p(k, t)$ that the system is in state $k$ at time $t$. The Chapman-Kolmogorov equations for the open model are then

$$ (1) \frac{d}{dt} p(k, t) = - (\lambda(K) + \sum_{i=1}^{n} \mu_i(k_i) (1-p_{ii})) p(k, t) \\ + \sum_{i=1}^{n} \lambda(K-1) p_{oi} p(b(k, i), t) \\ + \sum_{i=1}^{n} \mu_i(k_i+1) p_{i, n+1} p(a(k, i), t) \\ + \sum_{i \neq j}^{n} \sum_{i=1}^{n} \mu_i(k_i-1) p_{ij} p(C(k, i, j), t) $$

where it has been assumed that $p(k, t) = 0$ for any $k$ containing a negative element. We shall assume that $0 < \mu_i(k_i) < 1$ for all $i$, except for $k_i = 0$ where $\mu_i(0) = 0$ and that the Markov chain $\{\gamma_m\}$ has a single absorbing state $(n+1)$. Under these assumptions it can be shown that (1) represents the state transitions of an irreducible continuous time Markov chain. Foster's theorems [20]* guarantee the existence of a joint equilibrium probability distribution for the number of customers in each queue if there exists any positive solution $\{p(k)\}$ to (1) with $(dp(k, t)/dt) = 0$ such that $\sum p(k)= 1$.

<footnote>* Although Cohen's book [20] is not the primary reference for Foster's theorems the account given of them is very good.</footnote>
Let $\hat{P}$ be obtained from $P$ by deleting its first row and last column, and $e = (e_1, \dots, e_n)$ be a row vector such that

$$ (2) \qquad e = (p_{01}, \dots, p_{0n}) + e \hat{P} $$

Denote

$$ N(k) = \sum_{i=1}^{n} \sum_{m=1}^{k_i} \mu_i(m) $$

and

$$ \omega(K) = \sum_{m=0}^{K-1} \lambda(m) $$

$$ T(K) = \sum_{\substack{\text{all } k \\ \text{s.t. } \sum_{i=1}^{n} k_i = K}} N(k) $$

For the open model we have described, Jackson's theorem becomes

**Theorem 1** : The equilibrium joint probability distribution {$\hat{p}(k)$} for the number of customers exists if $e$ of (2) has a unique non-negative solution and $\sum_{K'=0}^{\infty} \omega(K')$ $T(K') < \infty$. It is given by

$$ p(k) = \left[ \sum_{K'=0}^{\infty} \omega(K') T(K') \right]^{-1} N(k) \omega(K) $$
The elements of e have a simple physical interpretation ; $e_i$ is the average number of times a customer transits through station i before leaving the system.

Now consider the "closed" model in which a constant number K of customers circulates within the QN indefinitely. The routing behaviour is given by a first-order irreducible Markov chain with transition matrix $P = (p_{ij})_{nxn}$. The Chapman-Kolmogorov equations are now

$$ (3) \quad \frac{d}{dt} p(k, t) = - \left( \sum_{i=1}^{n} \mu_i(k_i) (1-p_{ii}) \right) p(k, t) \\ \qquad + \sum_{i \neq j}^{n} \sum_{k_i=k}^{n} \mu_i(k_i+1) p_{ij} p(c(k, i, j), t) $$

where we have used the same notation as previously, including the rule that $p(k, t) = 0$ if k contains a negative element and that $0 < \mu_i(k_i) < \infty$ except $\mu_i(0) = 0$.

In this case Jackson's theorem becomes identical to a result discovered much later by Gordon and Newell [21]. We shall refer to Theorems 1 and 2 as the JGN theorems.

**Theorem 2** : Let $\pi$ be the unique probability row vector satisfying $\pi P = \pi$ and $\pi_i$, $1 \le i \le n$, its i-th element. Then $p(k)$, the equilibrium probability distribution satisfying (3) with $\sum_{i=1}^{n} k_i = K$, $k_i \ge 0$, exists and is given by

$$ p(k) = v(k) [G(K)]^{-1} $$

where

$$ v(k) = \sum_{i=1}^{n} \pi_i \frac{k_i}{\pi_i - \sum_{j=1}^{n-1} \pi_j} \frac{1}{\mu_i(k_j)} $$

and $G(K)$ in a normalizing constant so that $\sum p(k) = 1$.
Proofs of Theorems 1 and 2 are quite simple since they can be carried out by direct substitution in the equilibrium equations obtained by setting the time derivatives of $p(k; t)$ to zero in (1) and (3).

The JGN models present several shortcomings in their applicability to computer system modelling. The main ones are :

- (i) inability to distinguish among classes of customers with distinct stochastic behaviour

- (ii) restriction on service time and interarrival time distributions which must be exponential

- (iii) restriction on probabilistic routing behaviour given by a first order Markov chain

- (iv) inability to handle certain service disciplines of interest such as round-robin (RR)

- (v) exclusion of QN in which service time parameters depend on the number, and properties, of customers in a subnetwork rather than in each queue (includes blocking for example).

The remaining parts of this paper are devoted to results which allow us to remove most of these restrictions. There is at this time, however, no one single model for which equilibrium probability distributions are known and which allows the relaxation of all of the constraints (i) - (v).

The theorems of Jackson and Gordon and Newell have played and still play a central role in research on the theory of queueing networks. We shall see that more general network solutions can be often obtained as equivalent JNG networks. These transformations also serve as a proof technique in establishing properties of more general QN. Also, the diffusion approximation results presented in [6,9] are partly based on postulating a structure for the equilibrium probabilities similar to that presented in this section.
### III. WORK RATE THEOREMS

Exact solutions for the equilibrium probabilities of queueing networks with first-come-first-served (FCFS) service discipline and general service times are unavailable at present. The diffusion models presented in Section V yield only approximate solutions. Some interesting general results can, however, be given for the equilibrium work-rates of the service stations. These are mainly due to Chang and Lavenberg [17] and Lavenberg, Traiger and Chang [18].

Unless stated otherwise, the QN examined in this section are closed with K customers and n service stations. A customer proceeds to queue j after terminating service at queue i with probability $p_{ij}$. The matrix $P = (p_{ij})_{nxn}$ is irreducible. Service times at station i are independent and identically distributed with mean $O < \mu_i^{-1} < \infty$ and general distribution function $F_i = F_i(t)$; station i contains $m_i$ service centers with $m_i \ge 1$, $1 \le i \le n$; the case $m_i \ge K$ is possible and corresponds to a station with no queueing delay. If $F_i(t)$ is an exponential distribution, this assumption will correspond to taking $\mu_i(k_i)$ of the form

$$ \mu_i(k_i) = \begin{cases} k_i \mu_i & \text{if } 0 \le k_i \le m_i \\ m_i \mu_i & \text{if } k_i > m_i \end{cases} $$

Service times are independent from one station to another and are mutually independent at any given station; the service discipline is FCFS.

The state of the model is characterised by the number of customers and by the service time which has been obtained by those in service at each station. We define the service rendered by station i in $[0, t]$ as $W_i(t)$. The following theorem, due to Chang and Lavenberg [18], is stated without proof.
**Theorem 3** : $u_i = \lim_{t \to \infty} \frac{W_i(t)}{t}$ exists with probability one and is independent of the state at $t = 0$. It is the expected number of active servers station i at equilibrium.

The next result is a direct consequence of Theorem 2.

**Theorem 4** : If the model of this section is restricted to $F_i(t) = 1-\exp(-\mu_i t)$ then

$$u_i = \pi_i \mu_i^{-1} \frac{G(K-1)}{G(K)}$$

for $1 \le i \le n$.

**Proof** : We have

$$u_i = \sum_{k_i=1}^{m_i} k_i \frac{\Sigma_k p(k)}{k \Sigma_{k_i=m_i+1} p(k)} \frac{\Sigma_k p(k)}{k \Sigma_{k_i=m_i} p(k)}$$

yielding the result after some algebra.

**Theorem 5** : If $K = 1$ then for $1 \le i \le n$,

$$u_i = \frac{\pi_i \mu_i^{-1}}{\Sigma_j \pi_j \mu_j^{-1}}$$

This result is simply a direct consequence of semi-Markov process theory (see, for instance, Howard [22]). For $K=1$ the "position" of the customer is a stochastic process described by a first order Markov chain for state transitions and holding times in each state defined by general distribution functions. The work-rate of station i is then the equilibrium probability of finding the customer in that station.
$u_i$, which is the equilibrium value of the expected number of servers active at station i, can be interpreted as the equilibrium probability that station i is active when $m_i = 1$.

A result analogous to Theorem 4 has been established [17] for QN in which the only non-exponential service distribution are associated with stations at which no queueing can occur. This result can now be considered as a special case of the BCMP theorem [2] given in the next section.

**Theorem 6** : Maintaining the other assumptions stated at the beginning of this section suppose that $F_i(t) = 1-\exp(-\mu_i t)$, $1 \le i \le b$, and $m_i \ge K$ for $b+1 \le i \le n$. Then for all $1 \le i \le n$

$$u_i = \pi_i \mu_i^{-1} \frac{G(K-1)}{G(K)}$$

**Proof** : The idea of the proof is simple. We first show that Theorem 6 is true if each $F_j$, $b+1 \le j \le n$, is a mixture of Erlang distributions. Then we establish a well-known result : any distribution function of finite mean for a non-negative random variable can be approximated arbitrarily closely by a mixture of Erlang distributions. Under the assumptions of the theorem we are then able to invoke continuity of the work-rates to complete the proof. We shall proceed as in [17].

**Lemma 1** : If the servers at which no queueing can occur have a service time distribution which is a mixture of Erlang laws then Theorem 6 is true.

An Erlang-q distribution with mean $q/\alpha$ corresponds to the distribution function of a random variable which is the sum of q exponential random variables, each of parameter $\alpha$; it is therefore given by

$$E(t) = 1 - \sum_{j=0}^{q-1} \left( \frac{\alpha t}{j!} \right)^j e^{-\alpha t}$$
$$
\text{Now let } F_i(t) = \sum_{r=1}^{s_i} C_{ir} E_{ir}(t)
$$

where $E_{ir}(t)$ is the law $E(t)$ with parameters $q_{ir}$ and $\alpha_{ir}$, and $C_{ir} > 0$,
$$
\sum_{r=1}^{s_i} C_{ir} = 1
$$

We now construct a new QN by adding to <u>each</u> server of station i, $s_i$ new
stations i1,...,is_i succeeding the previous i-th station ; a customer
entering any server in station i in the old network will enter an exponen-
tial station i1 in the new network whose average service time is $\alpha_{ir}^{-1}$
with probability $C_{ir}$, $1 \le r \le s_i$. It will then proceed through $q_{ir}^{-1}$ other
identical stations before going to station j with probability $p_{ij}$. Notice
that the new QN has exponential service at all stations and no queueing at
stations derived from no-queueing stations of the old network. Each "new"
station sequence i1,..., is_i must contain only one customer : this is not
really a restriction since there are $m_i \ge K$ such sequences. Therefore
Theorem 2 is applicable. Call $k' = (k_1, ..., k_b, k_{b+1}, 1, ..., k_{b+1}s_1, ..., k_n s_n)$ and
P' the new transition matrix for routing behaviour. Let $\pi' P' = \pi'$ be the
new vector of equilibrium probabilities corresponding to P' : $\pi_i$ in $\pi$ has
been replaced by $\pi_{i1}, ..., \pi_{is_i}$ in $\pi'$ for $b+1 \le i \le n$. It is easily seen
that P' will also be irreducible, furthermore $\pi_i' = a\pi_i$, $1 \le i \le b$, and
$\pi_{ir}' = a C_{ir} \pi_i$ for $b+1 \le i \le n$ and $1 \le r \le s_i$. We then have,
$$
p(k') = [g(K)]^{-1} \sum_{i=1}^{b} \pi_i (a\pi_i \mu_i^{-1})^{k_i} \sum_{i=b+1}^{s_i} \pi_i \frac{(q_{ir} a C_{ir} \pi_i \alpha_{ir})^{-1}^{k_i}}{k_i!}
$$

However

$$
\begin{align*}
& \sum_{k_{i1}, \ldots, k_{is_i}, s.t. \Sigma k_{ir} = k_i} \\
& \qquad = \frac{1}{k_i!} \left( \sum_{r=1}^{s_i} z_{ir} \right)^{k_i}
\end{align*}
$$
and $\sum_{r=1}^{s_i} q_{ir} C_{ir} \alpha_{ir}^{-1} = \mu_i^{-1}$. Therefore

$$p(k) = a^k [f(k)]^{-1} v(k)$$

We normalize setting $[g(K)]^{-1} = a^{-K}[G(K)]^{-1}$ completing the proof since Theorem 6 is true if all service times are exponential.

The next lemma is classical (see Cox [24] for related matters).

**Lemma 2** : Let $F(t)$ be an arbitrary distribution function for a non-negative random variable of finite average. There exists for each $\varepsilon > 0$ a distribution function $\Phi(t)$ which is a convex combination of Erlang laws such that

$$\int_{0}^{\infty} |F(t) - \Phi(t)| dt < \varepsilon, \quad \left| \int_{0}^{\infty} t(F(t) - \Phi(t)) dt \right| < \varepsilon$$

**Proof** : We follow closely the proof in [17]. The idea is to construct a staircase approximation

$$\bar{F}(t) = \sum_{r=1}^{s} C_r u(t - rh)$$

to $F(t)$, where $h$ is the discretisation step, $u(t)$ is the unit step function and $s$ has been chosen so that

$$\int_{\tau}^{\infty} (1 - F(t)) dt < h$$

where $s$ is the integer part of $\tau/h$ plus 2. Also

$$C_r = F(h(r+1)) - F(h_r), \quad 2 \le r \le s-1$$

$$C_1 = F(2h), \quad C_s = 1 - F(sh)$$
Notice that $C_r \ge 0$, $\Sigma C_r = 1$. Each element step function is then shown to approximate an Erlang function establishing the correspondence between $F(t)$ and a convex combination of Erlangs. Clearly

$$ \int_0^\infty |F(t) - \bar{F}(t)| dt \le 2hC_1 + \sum_{r=2}^{s-1} hC_r + hC_s = 2h $$

Also

$$ \begin{aligned} & \int_0^\infty \left| 1 - \sum_{r=0}^{s-1} \left( \frac{\alpha t}{r!} \right)^r e^{-\alpha t} - u(t - s/\alpha) \right| dt \\ &= 2 \int_0^{s/\alpha} dt \left( 1 - \sum_{r=0}^{s} \left( \frac{\alpha t}{r!} \right)^r \right) = \frac{2se^{-r}r^r}{\alpha r!} \end{aligned} $$

which tends to zero for $r \to \infty$. Then

$$ \int_0^\infty |F(t) - \Phi(t)| dt \le \int_0^\infty |F(t) - \bar{F}(t)| + \int_0^\infty |\bar{F}(t) - \Phi(t)| < \varepsilon $$

The proof that the absolute value of the difference of the expected values is less than $\varepsilon$ is analogous.

To complete the proof of Theorem 6 an additional result is used from [18].

**Lemma 3** : $F = (F_1, \dots, F_n)$ be a set of probability distribution functions for a QN satisfying the assumptions of this section and let $\{F_i^m\}$ be a sequence of sets of such functions with

$$ \lim_{m \to \infty} \sum_{i=1}^{n} \int_{0}^{\infty} |F_{i}^{m}(t) - F_{i}(t)| dt = 0 $$
Then if $u_i(F)$ is the work-rate of station i with distribution functions $F_1, \dots, F_n$ we have

$$ \lim_{m \to \infty} u_i(F^m) = u_i(F) $$

From this result the theorem is now immediate. Notice that in the proof of theorem 6, and in particular in the proof of Lemma 1, our equations imply that stations 1 to b contain one server only ; clearly, the generalisation to $m_i > 1$ servers at these stations is immediate by modifying the expression for $p(k')$.

For the QN model this section Chang and Lavenberg have established the following :

**Theorem 7 [17]** : For the QN networks of this section

$$ \frac{u_i}{u_j} = \frac{\pi_i u_i^{-1}}{\pi_j u_j^{-1}} $$

This result can be of use in simulation studies of queueing networks, since it implies that work-rates for all queues can be deduced if the work-rate at only one station is known (or has been measured.

On the other hand, Lemma 3 can be helpful in computing work-rates approximately; if one replaces general service time distributions by appropriately chosen distributions which can be expressed as linear combinations of Erlang distributions, the resulting model is a continuous-time Markov chain which can (in principle) be solved numerically. Lemma 3 guarantees that in this manner one can compute work-rates with any desired degree of accuracy.

Another useful result has been established in [17].
**Theorem 8** : Let $u_i^{(K)}$ be the work-rate at station i, $1 \le i \le n$, of the closed QN with general service time distributions and K customers. Then

$$u_i^{(K)} \le u_i^{(K+1)}$$

for all $K = 0,1, ...$

This result has the simple consequence that Theorem 5 provides a lower bound to the values of the work-rates for each station, while upper bounds can be obtained by treating the case $K \to \infty$. In the limit $(K \to \infty)$ one of the stations (say j) will be saturated so that $u_j^{(\infty)} = 1$; all other $u_i^{(\infty)}$ can then be obtained by Theorem 7.
## XV → THE BCMP THEOREM

The major extension of the work of Jackson came a decade after his paper [1]. Although Jackson's work was in the line of traditional operations research, the BCMP theorem was motivated by problems of computer system modeling although it has far wider implications than just in computer science.

In this section we shall first consider the general model analyzed by Baskett, Chandy, Muntz and Palacios [2]. We then state their theorem.

Consider a network of queues with n service stations. Customers belong to one of ξ classes at any instant of time during their existence in the network. Routing in the network is described by a Markov chain with transition matrix $P = (p_{i,c;j,c'})$ where $0 \le i, j \le n+1, 1 \le c, c' \le \xi$; $p_{i,c;j,c'}$ is the probability that a customer of class c leaving station i will enter station j as a customer to class C' if $1 \le i, j \le n$. $p_{o,c;j,c'}$ is the probability that an external arrival first enter station j ($j \ne n+1$) in class c', and $p_{i,c;n+1,c'}$ is the probability that a customer of class c leaves the network after receiving service at station i.

The most general service time distribution considered is one which has rational Laplace transform; we shall call such service time distributions Coxian [24]. Notice that convex mixtures of Erlang distributions are Coxian so that, by Lemma 2 an arbitrary distribution may be approximated as closely as one may wish by a Coxian distribution. The service stations of the network are of one of the four following types.

**Type 1** : The station i has a single server with exponential service time of parameter $\mu_i(k_i)$ identical for all classes, $k_i$ being the total number

<footnote>+) A brief introduction to these distributions is given in the Appendix.</footnote>
of customers at the station, and the service discipline being FCFS.

**Type 2** : The service time distributions are Coxian and may be distinct for each class. There is a single server and the discipline is processor sharing (PS) so that each customer receives $k_i^{-1}$ seconds of service per second.

**Type 3** : The service time distributions are Coxian and may be distinct for each class. The number of servers at this station is larger than the maximum number of customers which may be queued here at any instant of time.

**Type 4** : The station contains a single server which allocates service according to a pre-emptive resume last-come-first-served discipline. That is, a new arrival interrupts the current service and takes over the server ; the interrupted customer resumes service when the customer which has interrupted it has finished . The interrupted customer's remaining service time is its total service time minus the amount it has already obtained. The service time distributions are Coxian and may be distinct for each class.

The arrival process to the QN is Poisson of parameter $\lambda(K)$ where K is the total number of customers in the network *)

In order to state the BCMP theorem we shall have to develop some notation. Let {$e_{jc}$} be a solution to the system of equations

$$e_{jc} = \sum_{c'=1}^{\xi} \sum_{i=1}^{n} e_{ic'} p_{i,c';j,c'} + p_{o,c';j,c}$$

for $1 \le j \le n$, $1 \le c \le \xi$. $e_{jc}$ is the relative frequency (or the normalized average number) of visits of a class c customer to service center j.

<footnote>*) A somewhat different arrival pattern is also allowed in [2].</footnote>
The state of the QN is the vector

$$ \hat{k} = (\hat{k}_1, \dots, \hat{k}_n) $$

where if station i is of

Type 1 :

$$ \hat{k}_i = (k_{i1}, \dots, k_{ik_i}) $$

$k_{ij}$ is the class of the j-th customer at station i in FCFS order.

Type 4 :

$$ \hat{k}_i = ( (k_{i1}, s_{i1}), \dots, (k_{ik_i}, s_{ik_i}) ) $$

$k_{ij}$ is the class of the j-th customer in LCFS order and $s_{ij}$ is its stage.

Type 2 or 3 :

$$ \hat{k}_i = ( (k_{i1}, s_{i1}), \dots, (k_{ik_i}, s_{ik_i}) ) $$

where $k_{ij}$, $s_{ij}$ are as in Type 4 except that the order is arbitrary. We denote by $k_{icm}$ the number of customers of class c in the m-th stage of service at station i.

We denote by $\mu_i(j)$ the parameter of the negative exponential service time distribution at a Type 1 station when $k_i=j$.

The Coxian densities will be represented as (see Figure 2) a sequence of exponential stages : $S_{ic}$ is the number of stages of the service time distribution of class c customers at station i and $\mu_{icm}$ is the parameter of the exponential distribution at its m-th stage. A customer's service time is then the sum of independent exponential random variables ; after the m-th stage the customer either heads for the $(m+1)$-th stage or ends its service with fixed branching probabilities. $A_{icm}$ is the probability that a class c customer at station i will reach its m-th stage of service.
For further details, it is suggested that the reader refer to the discussion
in the Appendix.

The BCMP Theorem is

**Theorem 9**: $p(\hat{k})$ the equilibrium probability distribution of $\hat{k}$ exists if $\gamma < \infty$. It is given by

$$p(\hat{k}) = \gamma^{-1} d(K) \prod_{i=1}^{n} f_i(\hat{k}_i)$$

where K is the total number of customers in state $\hat{k}$ and

$$f_i(\hat{k}_i) = \begin{cases} \frac{k_i}{\pi} \frac{e_{ik_i}}{\mu_i(j)}, & \text{if } i \text{ is Type 1} \\ j=1 & \\ \frac{\xi}{\pi} \frac{S_{ic}}{\pi} \frac{e_{ic} A_{icm} \mu_{icm}^{-1}}{\mu_{icm}^{\frac{k_i}{k}}}, & \text{if } i \text{ is of Type 2} \\ c=1 & \\ \frac{\xi}{\pi} \frac{S_{ic}}{\pi} \frac{e_{ic} A_{icm} \mu_{icm}^{-1}}{\mu_{icm}^{\frac{k_i}{k}}}, & \text{if } i \text{ is of Type 3} \\ c=1 & \\ \frac{k_i}{\pi} \frac{e_{ik_i} A_{ik} \mu_{ik}^{-1}}{\mu_{ik} \mu_{ij} s_{ij}}, & \text{, if } i \text{ is of Type 4} \\ j=1 & \end{cases}$$

where

$$d(K) = \begin{cases} \frac{K-1}{\pi} \lambda(m) & \text{if the QN is open} \\ m=0 & \\ 1 & \text{if it is closed} \end{cases}$$

$$\gamma = \sum_{\hat{k}} h(K) \prod_{i=1}^{n} f_i(\hat{k}_i)$$

The solution presented in Theorem 9 is for the equilibrium state probabilities which contain more detailed information than is often required (e.g. the ordering of customer classes in queue) and more information than is meaningful (e.g. the stage of service of customers). The following two corollaries concern marginal distributions obtained by aggregating states.
### Corollary 9.1

Consider a queueing network satisfying the conditions of Theorem 9. Let $\bar{k}_{ic}$ denote the number of customers of class c at station i. Let $\bar{k}_i = (\bar{k}_{i1}, \bar{k}_{i2}, \dots, \bar{k}_{i\xi})$ and $\bar{k} = (\bar{k}_1, \bar{k}_2, \dots, \bar{k}_n)$. Then the equilibrium probabilities of $\bar{k}$, $p(\bar{k})$, exist under the same conditions as $p(k)$ and are given by

$$p(\bar{k}) = \gamma^{-1} d(K) \prod_{i=1}^{n} g_i(\bar{k}_i)$$

where $\gamma$ and $d(K)$ are as in Theorem 9;

$$g_i(\bar{k}_i) = \begin{cases} \sum_{c=1}^{\bar{\kappa}_i} \frac{1}{k_{ic}!} e_{ic}^{k_{ic}} \left(\frac{1}{\mu_i}\right)^{k_i} & \text{if } i \text{ is of Type 1} \\ \sum_{c=1}^{\bar{\kappa}_i} \frac{1}{k_{ic}!} \left(\frac{e_{ic}}{\mu_i c}\right)^{k_{ic}} & \text{if } i \text{ is of Type 2 or 4} \\ \sum_{c=1}^{\bar{\kappa}_i} \frac{1}{k_{ic}!} \left(\frac{e_{ic}}{\mu_i c}\right)^{k_{ic}} & \text{if } i \text{ is of Type 3} \end{cases}$$

The most obvious implication of this corollary is that the distribution for $k$ depends only on the means of the service time distributions.

### Corollary 9.2

Consider an open network satisfying the conditions of Theorem 9 and with the additional constraints that the arrival process and service rates are constant (not state dependent). Let $k = (k_1, \dots, k_n)$. The equilibrium probabilities of $k$, $p(k)$, exist under the same condition as in corollary 9.1 and are given by

$$p(k) = \prod_{i=1}^{n} h_i(k_i)$$
where

$$h_i(k_i) = \begin{cases} (1-p_i)\rho_i^{k_i} & \text{if } i \text{ is of Type 1,2, or 4} \\ e^{-\rho_i} \frac{\rho_i^{k_i}}{k_i!} & \text{if } i \text{ is of Type 3} \end{cases}$$

and $\rho_i$ is the utilization of the i-th station given by $\sum_{c=1}^{\xi} e_{ic} \frac{1}{\mu_{ic}}$

The BCMP theorem contains as special cases stationary solutions to several models of interest. Consider, for instance, Jackson's QN in which routing is modified so as to be represented by a semi-Markov process [22].

We shall consider a QN identical to that of Section II except that, after completing its service at station i, a customer takes $\tau_{ij}$ seconds to arrive at queue j, where $\tau_{ij}$ is a random variable whose distribution function $F_{ij}(t)$ is Coxian. Of course, this customer will choose station j with fixed probability $p_{ij}$. This model is covered by the BCMP theorem if we add to the original QN $n^2$ new stations (i,j) of Type 3. The new routing probabilities will be

$$p'_{i}(i,j) = p_{ij}$$

$$p'(i,j),j = 1$$

and $p'_{o,i} = p_{o,i}$, $p'_{i,n+1} = p_{i,n+1}$ for all $1 \le i,j \le n$, considering an open QN. If $e'_i$ or $e'_j$ are the expected number of visits of a customer to a
station in the new QN and $e_i$ are the corresponding parameters in the original QN, we shall have $e'_{(i,j)} = e_i P_{ij}$, $e'_i = e_i$ for all $1 \le i,j \le n$. The equilibrium probability distribution of $k = (k_1, \dots, k_n)$ is then

$$p(k) = G \exp \left( \sum_{i=1}^{n} e_i \sum_{j=1}^{n} p_{ij} T_{ij} \prod_{m=0}^{K-1} \lambda(m) \prod_{i=1}^{n} e_i \prod_{j=1}^{k_i} \mu_i(j)^{-1} \right)$$

where $T_{ij} = \int_0^\infty t F_{ij}(t) dt$, and $k_i, 1 \le i \le n$, refers to the number in queue i. $G$ is the normalizing constant.

Models of this type, in which transitions of customers from one station (after completing their service) to another do not take place instantaneously can be of interest in representing certain aspects of multiprogramming systems. For instance, suppose a job terminates a period of service at the CPU because of an input-output request. Before this request can be placed in the queue of some input-output device it may be necessary to await the availability of input-output buffer space in primary memory. Similarly, before a certain job which has just terminated input-output can join the queue of "ready" jobs at the CPU it may have to wait until it can regain primary memory space (e.g. to load its "working-set" of pages). These can be considered as being random delays of the type discussed above.

Another apparent generalization of the BCMP theorem which is already contained in this important result concerns the generalization to routing probabilities defined by an h-order (rather than first order) Markov chain. Consider, for instance a job representing a user editing a text at a time-sharing terminal. The behaviour of the job could be characterized by a time spent at the terminal (thinking and typing) followed by a short execution at the CPU (interpreting a file request command), followed by
an input-output (loading a file into primary memory), followed by a CPU execution (modifying the contents of the file), and finally a return to the terminal. Clearly, the probability of transition of the job from one station to the other in this example can depend on, at least, the last two stations it has visited and not just the last one.

To see that this generalization can be reduced to the cases covered by the BCMP theorem, suppose that the routing probabilities of customers are given by the matrix $(p_{I,j;C,k})$ where $I = (i_1, i_2, ..., i_h)$ represents the last h stations visited by the customer $(i_h)$ is the most recent one) and j is the next station it will visit ; c and $c'$ are, as usual, the present and next class to which it will belong. We now obtain a new QN, equivalent to this one, but of the form covered by the BCMP theorem. For each class c and (h-1)-tuple of service stations $(i_1, ..., i_{h-1})$ of the original QN we form a class in the new QN of the form $(i_1, ..., i_{h-1}, c)$. The new transition probabilities are now :

$$p'_{i_h,K;j,K'} = p_{I,c;j,c'}$$

where

$$K = (i_1, ..., i_{h-1}, c), K' = (i_2, ..., i_{h-1}, i_h, c')$$

so that routing probabilities are now defined by a first-order Markov chain. The service time distributions for a customer of some class K in the new network are identical to those of a customer of the corresponding class C in the original QN. Otherwise the two networks are identical. Equilibrium probabilities, and other statistics, for the original network can be found by applying the BCMP theorem to the new network and summing over all values of $i_1, ..., i_{h-1}$ for a class K to obtain the corresponding quantity for class C.
Another useful special case of the BCMP theorem is obtained by considering closed subchains of the transition probability matrix $\underline{P} = (p_{i,c;j,c})$ for closed networks when $\underline{P}$ is completely decomposable. We may then write

$$ \underline{P} = \begin{bmatrix} \underline{P}_1 & & \\ & \underline{P}_2 & 0 \\ & & \ddots \\ & & & \underline{P}_L \end{bmatrix} $$

where $\underline{P}_\ell$ is stochastic and non-decomposable and represents an ergodic Markov chain with state space $E_\ell$, $1 \le \ell \le L$. Since there is no exchange of customers between each of the sub-chains $\underline{P}_1, \dots, \underline{P}_L$, the number of customers associated with each sub-chain remains constant in a closed network. The following result shows that at steady state the marginal probability distributions associated with this model are identical to the probability distributions of a model with the same number of stations and L customer classes where the parameters of each customer classes are obtained from parameters related to the closed sub-chains, and in which no class changes occur.

Theorem 11 (Marginal distributions with closed sub-chains) Let $k_i'$ represent the total number of customers at station i belonging to the $\ell$th non-decomposable sub-chain $\underline{P}_\ell$ and

$$ \hat{k}'_i = (k'_{i1}, \dots, k'_{iL}) $$

Consider the stationary joint probability distribution $p(\hat{k}')$ of finding station i in state $\hat{k}'_i$, $1 \le i \le n$, where

$$ \hat{k}' = (k'_1, \dots, k'_n) $$

Denote

$$ \rho_{il} = \sum_{c \text{ s.t. } (i,c) \in E_\ell} e_{ic} \mu_{ic} $$

$$ e_{il} = \sum_{c \text{ s.t. } (i,c) \in E_\ell} e_{ic} $$
if i is of Type 2,3,4. Then for the closed queueing network described above

$$p(\hat{k}') = \alpha^{-1} \prod_{1}^{n} f_{i}(\hat{k}'_{i})$$

where α is the normalizing constant and

$$f_i(\hat{k}'_i) = \begin{cases} \prod_{j=1}^{k_i} \mu_i(j)^{-1} \prod_{\ell=1}^{L} \frac{1}{k_{i\ell}'} \hat{e}_{i\ell}^{k_{i\ell}'}, & \text{if } i \text{ is of Type 1} \\ \prod_{\ell=1}^{L} \frac{1}{k_{i\ell}'} \rho_{i\ell}^{k_{i\ell}'}, & \text{if } i \text{ is of Type 3} \\ \prod_{\ell=1}^{L} \frac{1}{k_{i\ell}'} \rho_{i\ell}^{k_{i\ell}'}, & \text{if } i \text{ is of Type 2 or 4.} \end{cases}$$

Theorem 11 allows the computation of stationary probabilities of the number of customers of each sub-chain at the service stations; the form of the solution is identical to the case where each sub-chain is viewed as a customer class with appropriately chosen parameters.

The four types of service centers included in the BCMP theorem are not all motivated by practical considerations. However Type 2 service centers are, in many cases, a reasonable representation of central processing units allocated in quanta; the processor sharing discipline is an idealization with a quantum of "very small" (in fact, zero) duration and no overhead associated with switching from one job to the other. A Type 3 station represents well terminals in a time-sharing system. Type 1 service centers are useful in many instances (secondary memory units, input-output devices, etc.). There are no concrete examples for the use of Type 4 stations in computer systems.
V. POISSON PROCESSES AND OUEUEING NETWORKS

The models discussed in the previous sections do not allow analytical
solution for certain apparently simple systems such as networks of queues with
first-come-first-served discipline and general service time distributions. Also,
even though an abundant literature exists on the subject, queues with generally
distributed and independent interarrival times and general service time
distributions are not completly amenable to analytical treatment. Fortunately,
Poisson processes may tend to appear in queueing networks due to sampling
or superposition, or because some queues transform Poisson arrival processes to
Poisson departure processes. We present here certain results of this nature
since they can serve as guidelines in order to make reasonable approximations
concerning Poisson arrivals to queues, and also to clarify some of the properties
of the networks of queues covered by the BCMP theorem.

The first result we established is easy to show though we have not been
able to find it in the classical literature.

Theorem 12 (Infrequent samplings of a renewal process).

Let $T = \{t_i \mid t_o = 0, t_{i+1} > t_i, i \ge 0\}$ be a renewal process; i.e.

$t_{i+1} - t_i$ are independent and identically distributed with common density function $f(t)$ (one may allow $t_1$ to have a different density function $f_1(t)$).

For each $t_i$ let $X(t_i)$ be a binary random variable such that $P_r[X(t_i) = 1] = p$
for all $i$; also $X(t_i)$ is independent of $X(t_j)$ for $i \neq j$. Consider the process
$\hat{T} = \hat{t}_i \hat{t}_{i+1} > \hat{t}_i \hat{t}_{i+1}$, $i>0$ where $\hat{t}_i$ is the i-th element of the sequence $T$ for
which the corresponding value of $X(.)$ is one. Let $m_j = \int^\infty \hat{t}_j^j f(t) dt$, $j=1,2,...$.
Under the following limiting conditions $\hat{T}$ tends to a Poisson process of parameter
$\lambda \hat{p}$ ($\lambda=\hat{m}_1^{-1}$):

(i) $\lambda \to \infty$; $p \to 0$ so that $\lambda p$ tends to a finite constant

(ii) $\lambda m_j \to 0$ for all $j>1$.

Proof. It is easily seen that $\hat{T}$ is a renewal process. Let $N_t$, $\hat{N}_t$ be the number of renewals up to time $t>0$ for $T$ and $\hat{T}$, respectively, with

$$
G(t,x) = \sum_{r=0}^{\infty} x^r \mathrm{Pr}[N_t=r]
$$

$$
\hat{G}(t,x) = \sum_{r=0}^{\infty} x^r \mathrm{Pr}[\hat{N}_t = r]
$$
and corresponding Laplace transforms $G^*(s,x)$, $\hat{G}^*(s,x)$. Then it is well known that

$$ \hat{G}^*(s,x) = \frac{1-f^*(s)}{s[1-xf^*(s)]} $$

if $f_1(t) = f(t)$, where $f^*(s)$ is the Laplace transform of $f(t)$. $f_1^*(s)$, the Laplace transform of the renewal density for $T$, is

$$ \begin{aligned} f_1^*(s) &= \sum_{n=1}^{\infty} (1-p)^{n-1} p[f^*(s)]^n \\ &= \frac{pf^*(s)}{1-(1-p)f^*(s)} \end{aligned} $$

and

$$ G^*(s,x) = \frac{1 - f^*(s)}{s[1-(1-p)f^*(s)-xpf^*(s)]} $$

Since the function $f^*(s)$ is analytic in the half-plane $Re(s)>s_0$, for some $s_0<0$, we can expand it about $s=0$:

$$ f^*(s) = 1 - m_1 s + m_2 s^2 / 2! - \dots $$

Therefore

$$ G^*(s,x) = \frac{m_1 s - m_2 s^2 / 2! + \dots}{s[m_1 s - m_2 s^2 / 2! + \dots + pf^*(s)(1-x)]} $$

Dividing both numerator and denominator by $m_1 s$ gives

$$ G^*(s,x) = \frac{1 - \lambda m_2 s^2 / 2! + \dots}{[s - \lambda m_2 s^2 / 2! + \dots + \lambda p f^*(s)(1-x)]} $$

Taking limits as in (i) and (ii) we remain with

$$ G^*(s,x) \rightarrow \frac{1}{s + \lambda p(1-x)} $$

But the right-hand term in the above expression is the value of $\hat{G}^*(s,x)$ with $f^*(s) = \lambda p/(\lambda p+s)$ i.e. for a Poisson process; this completes the proof.

One can think of many applications of this result to computer systems. For instance, suppose that page transfer requests to a paging drum arrive according to a renewal process; no available theoretical models allow us to treat this case, exactly.
If the number of sectors is large, however, and if requests are directed to the sectors with fixed probability equal for each sector we would expect that arrivals to each sector would be Poisson so that the system can be analyzed.

The next result is classical; it explains many observations of Poisson processes (or close approximations of Poisson processes) in telephone systems, computer networks, computer systems, etc.

**Theorem of PALM-KHINTCHINE [20].** Let $T^1, \dots, T^n$ be n independent renewal processes:

$$T^i = \left\{ t_j^i \mid j=0,1,\dots,t_0=0, t_{j+1}^i > t_j^i \right\}, \quad 1 \le i \le n$$

of rate $\lambda_1, \dots, \lambda_n$ respectively. Let $T$ be the renewal process such that $t_j$ is an element of $T$ if and only if it is an element of one of the $T^i$, $1 \le i \le n$.

Let $\lambda = \sum_{1}^{n} \lambda_{i}$. The process $T^*$ obtained from $T$ by letting $n \to \infty$, $\lambda_i \to 0$ for each $1 \le i \le n$ and $\lambda \to \lambda^*$ (a finite constant) is Poisson of rate $\lambda^*$. (The *rate* of a renewal process is the inverse of the expected time between successive renewals).

It is easy to show that in the limit $(t_{j+1}-t_j)$ is exponentially distributed with the same distribution for all $j>0$; it is more difficult to show that $T^*$ is a renewal process. The first part of the proof was established by PALM and the second part (the complete proof) by KHINTCHINE (see [20]); generalizations of this result continue to appear in the recent literature, in particular to cases where $T^1, \dots, T^n$ are not mutually independent.

This result indicates that we should expect to observe Poisson arrival processes to resources in computer systems which receive requests from a large number of devices operating independently of each other.

Now suppose that a queueing system receives requests for service according to a Poisson process; it is obviously interesting to determine under what conditions the departure process is Poisson. If this property holds we shall say that the system is *complete* (this is the $M \to M$ property discussed in [3]). Let $s(t)$ be the state at time $t$ of a Markovian queueing system; e.g. for the $M/G/1$ queue $s(t) = (N(t), X(t))$ where $N(t)$ is the number of customers in the system and $X(t)$ is the total service time already obtained by the customer in service at time $t$. Assuming a Poisson arrival process of rate $\lambda$ to the queueing system, we say that it is *complete* if as $t \to \infty$

$$\forall s, \Pr[i \text{ departures in } (t-\Delta t, t) \mid s(t) = s] = \begin{cases} \lambda \Delta t + o(\Delta t), & \text{if } i=1 \\ o(\Delta t), & \text{if } i \ge 1 \\ 1-\lambda \Delta t + o(\Delta t) & \text{if } i=0 \end{cases}$$
Let us denote

$$ \alpha(s',s,t) = \lim_{\Delta t \to 0} P_r [s(t+\Delta t)=s \mid s(t) = s'] / \Delta t $$

$$ \beta(s,s', t) = \lim_{\Delta t \to 0} P_r [s(t-\Delta t)=s' \mid s(t) = s] / \Delta t $$

If $p(s,t)$ is the probability of finding the system in state $s$ at time $t$, then clearly

$$ \beta(s,s',t) = \alpha(s',s,t) p(s',t) / p(s,t) $$

Now let $A(s)$ be the set of all states such that a state transition from $s' \in A(s)$ to $s$ is accompanied by a departure from the queueing system. The definition given above for a complete system can now be written as

$$ \forall s, \lim_{t \to \infty} \sum_{s' \in A(s)} \beta(s, s', t) = \lambda \quad \therefore \lambda = 1 $$

or

$$ \forall s, \lim_{t \to \infty} \sum_{s' \in A(s)} \alpha(s',s,t) p(s',t) = \lambda p(s,t) $$

Assuming that state transition rates are independent of time $[\alpha(s',s,t) = \alpha(s, s)]$ and that the equilibrium state probabilities exist $[\lim_{t \to \infty} p(s,t) = p(s)]$ then we have the condition

$$ \forall s \quad \sum_{s' \in A(s)} \alpha(s',s)p(s') = \lambda p(s) $$

When there are different classes of customers the queueing system is said to be complete if the departure process for each class of customers is Poisson when the arrival process for each class of customers is Poisson. Let $A_r(s)$ denote the set of all states such that a state transition from $s' \in A_r(s)$ to $s$ is accompanied by a departure of a class $r$ customer from the queueing system. The condition for completeness becomes

$$ (\forall s) \quad (\forall \text{customer class } r) \quad \sum_{s' \in A_r(s)} \alpha(s',s)p(s') = \lambda_r p(s) $$

**Theorem 13 (Completeness of BCMP systems)** The four types of service stations considered in the BCMP theorem are complete.

The BCMP theorem gives a solution for a class of queueing networks with four types of service stations. The solution has the so called "product form" i.e. there is one factor for each service center. It can be shown [3] that a
sufficient condition for a product form solution is that the service centers
all be complete.

**Theorem 14 (Product form solution and completeness) [3]**

Let N be a network of queues model in which each service station is
complete. Let ξ be the number classes of customers and n the number of service
stations. Let e<sub>ic</sub> be the mean arrival rate of class c customers to service
center i. Let {P<sub>i</sub>(S<sub>i</sub>)} be the equilibrium state probabilities for service center
i when the arrival process for class c customers is Poisson with mean rate
e<sub>ic</sub>, 1<c<ξ.

(a) If N is a closed network then the equilibrium state probabilities for the network are given by

$$
P(S_1, S_2, \dots, S_N) = \gamma^{-1} P_1(S_1) P_2(S_2) \dots P_N(S_N)
$$

where $\gamma^{-1}$ is a normalizing constant chosen so that the network state
probabilities sum to one.

(b) If N is an open network then the equilibrium state probabilities for the network are given by

$$
P(S_1=S_2=\dots=S_N) = P_1(S_1) P_2(S_2) \dots P_N(S_N)
$$

While this only a sufficient condition there are no solutions known to
the authors for networks of queues with general transition probabilities and
service stations which are not complete. (Another method of approaching the
search for solutions is "local balance" introduced by Chandy [35]. The relation-
ship between local balance, Poisson departures and the product form solution
is further explored by Chandy, et. al. in [37].)
## VI. Computational Algorithms for Closed Queueing Networks

The solutions for the equilibrium state probabilities for a large class of queueing network models have been presented in previous sections. These solutions are useful only if they allow the computation of various performance measures (utilizations, queue length distributions, etc.) at a computational cost which is competitive with that of simulation. For open network models the solutions involve infinite sums and therefore can be evaluated exactly only when the infinite sums have closed form solutions. For closed network models the amount of computation is always finite but unless some care is taken the computation can quickly become impractical as the number of customers and number of service stations increases. For example the normalization constant T(K) of Theorem 1 has been defined as a sum of terms. The number of terms in this sum is $\binom{K+n-1}{n-1}$ which, even for relatively small values of n and K, can be impractically large. Fortunately, methods have been developed for computing quantities of interest which are computationally much more efficient than the straightforward, brute-force approach. Most of these methods are originally due to Buzen [23]. Here we follow the "generating function approach" given by Williams and Bhandiwad [36] because it provides an orderly and unified approach to the problem. To illustrate this method we will develop the computational algorithms for a case of intermediate complexity: a closed queueing network with a single class of customers, single server service stations and state dependent service rates. (More details and the treatment of more general cases can be found in [36].)

For the class of network models described above the equilibrium state probabilities have the form:

$$P(k) = \frac{1}{C} \prod_{i=1}^{n} A_{i}(k_{i})$$
where
$$A_{i}(k_{i}) = \prod_{j=1}^{k_{i}} x_{ij}$$
$$x_{ij} = \lambda_{i}(\mu_{i}(j))^{-1}$$

where the normalizing constant C is

$$C = \sum_{\sum k_i = K} \prod_{i=1}^{n} A_i(k_i)$$
$$k_i \ge 0$$
**Notation:** The following notational conventions are used throughout this section.

A number of _generating functions_ are defined which are polynomials in a dummy variable $t$. Small letters will be used (perhaps with subscripts) to denote these functions (e.g. $g(t)$, $g_i(t)$). Capital letters are used to denote the coefficients in the generating functions. For example: $g(t) = G(0)+G(1)t+G(2)t^2...$

Now define the generating function for service station i, $a_i(t)$ as

$$a_i(t) = 1 + A_i(1)t + A_i(2)t^2 + A_i(3)t^3 + ...$$

The generating function for the network model is defined as

$$g(t) = \prod_{i=1}^{n} a_i(t)$$

The normalization constant $C$ is easily seen to be the coefficient of $t^K$ in $g(t)$. Since it is the simplest application we first show how $C$ may be calculated.

First define:

$$ (4) \qquad \begin{aligned} g_1(t) &= a_1(t) \\ g_i(t) &= \prod_{j=1}^{i} a_j(t) \end{aligned} $$

Clearly $g_n(t) = g(t)$.

Now, given any two generating functions

$$r(t) = R(0)+R(1)t+R(2)t^2+...$$

$$\text{and } q(t) = Q(0)+Q(1)t+Q(2)t^2+...$$

we define

$$R(j)*Q(0) = \sum_{\ell=0}^{j} R(j-\ell)Q(\ell)$$

Clearly $R(j)*Q(0)$ is the coefficient of $t^j$ in $r(t)q(t)$.
From Eq. (4) it is clear that

$$ (5) \qquad \begin{aligned} G_l(j) &= A_l(j) \\ G_i(j) &= G_{i-1}(j) * A_i(0) \quad , \quad 1 \le i \le n \end{aligned} $$

To calculate $C= G(K) (=G_n(K))$ one uses (5) to calculate the $G_l(j)$, $G_2(j)$, $1 \le j \le K$, etc. until $G_n(K)$ is reached. It is not difficult to show that this computation grows as $nK^2$.

Now consider the problem of computing the utilizations (i.e. steady-state probability that the server is busy) for the various servers. First define:

$$ h_i(t) = \frac{g(t)}{a_i(t)} (a_i(t)-1) $$

i.e. $h_i(t)$ is $g(t)$ with $a_i(t)$ replaced by $(a_i(t)-1)$. Thus the $t^0$ term is removed from the generating function $a_i(t)$.

Note that the coefficient of $t^K$ in $h_i(t)$ is given by

$$ H_i(K) = \sum_{\substack{n \\ \sum k_j = K}} \prod_{j=1}^{n} A_j(k_j) $$

Thus the utilization of the i-th server is just

$$ \frac{H_i(K)}{G(K)} $$

Now consider the calculation of the coefficients $H_i(K)$.

We can write

$$ \begin{aligned} h_i(t) &= g(t)[1-d_i(t)] \\ &= g(t)-g(t)d_i(t) \end{aligned} $$

where $d_i(t)$ is the inverse of $a_i(t)$, i.e. $a_i(t)d_i(t) = 1$

Clearly then

$$ H_i(j) = G(j) - [G(j) * D_i(0)] \quad (6) $$

Thus, once the coefficients $D_i(t)$ are known, this computation is straightforward and requires on the order of $nK^2$ steps to compute $H_i(K)$, $1 \le i \le n$.
The $D_i(l)$ are easily determined from the condition $d_i(t)a_i(t) = 1$ [36].

One obtains:

$$D_i(0) = 1$$ $$D_i(l) = \sum_{v=1}^{l} D(l-v)A(v) \quad (7)$$

This computation grows as $K^2$.

Summarizing the above we have the following procedure for calculating utilizations.

1. Compute $D_i(l)$ $1 \le i \le n, 0 \le l \le K$ from Eq. 7

2. Compute $H_i(K)$ $1 \le i \le n$ from Eq. 6

3. Compute $G(K)$ from Eq. 5

Since the computation of the $D_i(l)$ grows as $K^2$, the computation of the $H_i(K)$ as $nK$ and the computation of $G(n)$ as $nK^2$, the computation of all utilizations grows as $nK^2$.

In special cases the amount of computation can be reduced considerably. For example, take the case in which the service rates are not state dependent

$(\mu_i(j) = \mu_i, \forall j>0)$. Then

$$A_i(k_i) = X_i \quad \text{where } X_i = \frac{\lambda_i}{\mu_i}.$$

In this case $a_i(t) = 1+X_i t+X_i^2 t^2+...$

Therefore

$$g_i(t) = g_{i-1}(t)a_i(t) = g_{i-1}(t)\left[\frac{1}{1-X_i t}\right]$$

$$g_i(t) = g_{i-1}(t) + X_i \operatorname{tg} i(t)$$

Clearly then:

$$G_i(j) = G_{i-1}(j) - X_i G_{i-1}(j-1)$$

Using this expression for calculating $G(K)$ the computation grows as $nK$ rather than $nK^2$.
The $H_i(K)$ are more efficiently calculated in this case as well. As defined previously

$$h_i(t) = g(t)[1-d_i(t)].$$

But now $d_i(t)$ is clearly $1-X_i t$ so that we can write

$$h_i(t) = X_i \operatorname{tg}(t)$$

$$\text{Clearly then } H_i(n) = X_i G(n-1)$$

and the utilization of the i-th server is simply

$$\frac{X_i G(n-1)}{G(n)}$$
# VII- CONCLUSIONS

In this paper we have presented a survey of exact results concerning networks of queues which are applicable to the evaluation of the performance of complex computer systems.

The results summarized here include expressions for equilibrium probability distributions for queues with exponential service time distributions or with service time distributions having rational Laplace transforms and specific service disciplines. The routing behaviour of the customers in these networks can be described by arbitrary (but finite) order Markov chains or by semi-Markov processes, and classes of customers possessing distinct statistical behaviour are allowed.

For networks of queues for which equilibrium probability distributions are unavailable, we have presented work-rate theorems relating the utilisation factors of the different service stations. The networks for which these results are given (general service time distributions at the different service stations and routing described by a first order Markov chain) have every chance of not being solved exactly in the future, because special instances have been examined with very limited success since the inception of queueing theory. Therefore in any application, such models will have to be analyzed via a Monte-Carlo simulation; the work-rate theorems are important in this respect since they allow the deduction of all work-rates from measurements of only one of them. They can be also useful in choosing best estimators for a given quantity. Several interesting problems of mathematical nature remain to be solved in this area.

Generalizations of the work-rate theorems to several classes of customers and to arbitrary order Markov chain representations of transition probabilities would be useful.
Also one might ask whether analogous theorems for equilibrium proba-
bilities might be established.

The most interesting direction for further research lies probably
in the study of approximate (e.g. diffusion or iterative) methods for solving gene-
ral classes of models. Finally studies related to the optimal and adaptative
control of models of multiprogramming systems [29, 38] will provide both simulatiny
theoretical work and practical applications to the design of self-regulating
operating systems.
# APPENDIX

The purpose of this section is to provide a brief introduction to Coxian or almost general probability distributions, that is those having rational Laplace transform. The presentation given here stresses the intuitive aspects.

Consider the network presented on Figure 4, and resume that between points A and B a single customer is present at a time; the time spent by the customer in this network will correspond to the realization of a random variable whose distribution is Coxian. The customer entering at point A either immediately leaves the system from B with probability ($1 - b_o$) or enters into a first stage with probability $b_o$, where it will remain for a time which is a realization of a random variable $\tau_1$ exponentially distributed with parameter $\lambda_1$; In general, after residing in the i-th stage ($1 \le i \le n-1$) the customer will either proceed to the (i + 1) - th stage with probability $b_i$ or leave through B with probability ($1 - b_i$); the time spent in stage i by the customer is exponentially distributed with parameter $\lambda_i$.

Let $f(t)$ be the probability density function for the time $\tau$ spent by a customer in traversing the network of Figure 4 from A to B and let $f^*(s)$ be its Laplace-Stieltjes transform :

$$f^*(s) = \int_{0}^{\infty} e^{-st} f(t) dt$$

Nothing that with probability ($1 - b_o$), the random variable $\tau$ has value zero, and that with probability $b_o$ it is the sum of a random number of exponential random variables, we have

$$f^*(s) = (1 - b_o) + \sum_{i=1}^{n} (1 - b_i) a_i \prod_{j=1}^{i} \frac{\lambda_j}{s + \lambda_j}$$

where

$$a_i = b_0 b_1 \cdots b_{i-1}, \quad 1 \le i \le n$$
Also

$$E\{\tau\} = \sum_{i=1}^{n} \frac{a_i}{\lambda_i}$$

and the variance of τ is

$$Var(\tau) = \sum_{i=1}^{n} \frac{a_i}{\lambda_i^2}$$

The class of density functions whose Laplace-Stietjes transform is of this form contains both hypoexponential and hyperexponential densities; in fact it contains all densities whose transform is a rational function of s and is closed under addition. Therefore it contains, in particular, the linear combinations of Erlang densities which have been used in the approximation results concerning work-rates presented in Section III.
<figure><img src="image_3.png" /><figcaption>Figure 1</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2</figcaption></figure>
REFERENCES

1. J.R. JACKSON, "Jobshop-Like Queueing Systems", Management Science 10, 131-142, 1963.

2. F. BASKETT, M. CHANDY, R. MUNTZ, J. PALACIOS, "Open, Closed and Mixed Networks of Queues with Different Classes of Customers", JACM 22, 2 April 1975.

3. R.R. MUNTZ, "Poisson Departure Processes and Queueing Networks", IBM Research Report RC4145, December 1972.

4. D.P. GAVER, "Probability Models for Multiprogramming Computer Systems", JACM 14, 423-438, 1967.

5. D.P. GAVER, G.S. SHEDLER, "Processor Utilization in Multiprogramming Systems Via Diffusion Approximations", Operations Research 21, 569-576, 1963.

6. H. KOBAYASHI, "Application of the Diffusion Approximation to Queueing Networks: Part I - Equilibrium Queue Distribution", JACM 21, 2, 316-328, April 1974.

7. E. GELENBE, "On Approximate Computer System Models", JACM 22, 2, 261-269, April 1975.

8. E. GELENBE, "A Non-Markovian Diffusion Model and its Application to the Approximation of Queueing and Computer System behaviour", Research Report 74-1, Chaire de Systèmes Informatiques, University of Liege, Belgium, 1974.

9. M. REISER, H. KOBAYASHI, "Accuracy of a Diffusion Approximation for Some Queueing Systems", IBM J. Res. Develop. 18, 2, 110-124, March 1974.

10. M. BADEL, A.V.Y. SHUM; "Accuracy of a Diffusion Approximation Method for Computing Work-Rates in Multiprogramming Systems", IRIA Research Report, 1975.

11. M. CHANDY, U. HERZOG, L. WOO, "Approximate Analysis of General Queueing Networks", IBM Research Report RC4931, July 1974.

12. P.J. COURTOIS, "On the Near Complete Decomposability of Networks of Queues and of Stochastique Models of Multiprogramming Systems", Carnegie-Mellon University, Computer Science Department, Research Report, November 1971.

13. P.J. COURTOIS, P. GEORGES, "An Evaluation of the Stationary Behaviour of computations in Multiprogramming Computer Systems", Proc. ACM-ICS 70, Vol. 1, 98-115, Bonn, 1970.

14. P.J. COURTOIS, "Decomposability, Instabilities, and Saturation in Multiprogramming Systems", Comm. ACM 18, 7, 371-377, July 1975.

15. B. AVI-ITZHAK, P. HEYMAN, "Approximate Models of Multiprogramming Computer Systems", Operations Research, 21, 1212-1230, 1973.

16. M. REISER, H. KOBAYASHI, "Blocking in a Queueing Network with Two Exponential Servers", IBM Research Report RJ1360, 1974.
17. A. CHANG, S.S. LAVENBERG, "Work-rates in Closed Queueing Networks with General Independent Servers", IBM Research Report RJ989, 1972.

18. S.S. LAVENBERG, I.L. TRAIGER, A. CHANG, "Bounds on Work-Rates for Closed Queueing Networks with Non-Exponential Stages", IBM Research Report RJ1182, March 1973.

19. W. FELLER, "Diffusion Processes in One Dimension", *Trans. Am. Math. Soc.*, 77, 1-31, 1954.

20. J. COHEN, "The Single Server Queue", North-Holland, Amsterdam, 1969.

21. W.J. GORDON, G.F. NEWELL, "Closed Queueing Systems with Exponential Servers", *Opns. Res.* 15, 254-265, 1967.

22. R.A. HOWARD, "Dynamic Probabilistic Systems, Vol. I and II", John Wiley, New York, 1973.

23. J. BUZEN, "Computational Algorithms for Closed Queueing Networks with Exponential Servers", *Comm. ACM* 16, 527-531, 1973.

24. D.R. COX, "A Use of Complex Probabilities in the Theory of Stochastic Processes", *Proc. Cambridge Philosophical Society* 51, 2, 313-319, April 1955.

25. D.P. GAVER, "Diffusion Approximations and Models for Certain Congestion Problems", *J. Appl. Prob.*, 5, 607-623, 1968.

26. G.F. NEWELL, "Applications of Queueing Theory", Chapman and Hall, Ltd., London, 1971.

27. J.P. BUZEN, "Queueing Network Models of Multiprogramming", Ph.D. Thesis, Harvard University, 1971.

28. E. GELENBE, "Modèles de Systèmes Informatiques", Thèse de Docotorat d'Etat, Université Paris VI, 1973.

29. A. BRANDWAJN, J. BUZEN, E. GELENBE, D. POTIER, "A Model of Performance for Virtual Memory Systems", Proc. ACM-SIGMETRICS Symposium, Montreal, September 1974.

30. R.R. MUNTZ, "Analytical Modelling of Interactive Systems", *Proceedings IEEE*, Special Issue on Interactive Systems, June 1975.

31. M. BADEL, E. GELENBE, J. LEROUDIER, D. POTIER, "Adaptive Optimization of a Virtual Memory System", *Proceedings IEEE*, Special Issue on Interactive Systems, June 1975.

32. B. COCHI, "Parameter Analysis of Some Queueing Models", Rapport de l'Ecole des Mines de Paris, February 1975.

33. G. FISHMAN, "Concepts and Methods in Discrete Event Digital Simulation", John Wiley, New York, 1973.
34. J. LEROUDIER, M. PARENT, "Quelques Aspects de la Modélisation des Systèmes Informatiques par Simulation à Evènements Discrets", to appear in *Revue d'Informatique de l'AFCET*.

35. K.M. CHANDY, "The Analysis and Solutions for General Queueing Networks", Proc. Sixth Annual Princeton Conf. on Information Sciences and Systems, Princeton, N.J., March 1972.

36. A.C. WILLIAMS, R.A. BHANDIWAD, "A Generating Function Approach to Queueing Network Analysis of Multiprogrammed Computers", to be published in *Networks*.

37. K.M. CHANDY, J. HOWARD, T.W. KELLER, D.J. TOWSLEY, "Local Balance, Robustness, Poisson Departures and the Product Form in Queueing Networks", submitted to *JACM*.

38. P.J. DENNING, S. GRAHAM, "Multiprogrammed memory management", Proc. *IEEE*, Vol. 63, N° 6, pp. 924-939, June 1975.