RIA
laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B.P. 5 78150 - Le Chesnay
France
Tél.: 954 90 20

F

laboratoire de recherche
en informatique
et automatique

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

A NON-MARKOVIAN
DIFFUSION MODEL
AND ITS APPLICATION
TO THE APPROXIMATION OF
QUEUEING SYSTEM BEHAVIOUR

Erol.

E. Gelenbe

Rapport de Recherche N° 158

Mars 1976
# A NON-MARKOVIAN DIFFUSION MODEL AND ITS APPLICATION TO THE APPROXIMATION OF QUEUEING SYSTEM BEHAVIOUR

E. Gelenbe

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

## Résumé :

Par l'inclusion de temps d'arrêt sur les frontières distribuées selon des lois générales, nous éten-
dons un modèle de diffusion au cas non-Markovien. Le modèle est appliqué au calcul des perfor-
mances de modèles de systèmes informatiques.

## Abstract :

We generalize a diffusion model of FELLER by the inclusion of holding times on the boundaries which are distributed according to a general probability distribution function. The stationary solution associated with the model is shown to depend only on the first moments of the holding time densities. The model is applied to the approximation of computer system behaviour.
# 1. Introduction

Models of multiple resource computer systems [2] have proved to be a good source of new problems and of general solution techniques [3-5] applicable to the analysis of networks of queues [1]. The research reported in this paper was motivated by this class of problems.

A promising method for the approximation of queueing systems with general service time distributions has originated with the work of GAVER [6] and NEWELL [7] who suggested the use of a diffusion process to approximate the number in queue. Consider for instance the GI/G/1 queue; basic to the diffusion approximation for this model is the assumption that as soon as a busy period begins (i.e. a customer arrives to a previously empty system) the stochastic process representing the number in queue is adequately approximated by the predictions of the central limit theorem which in reality are only valid asymptotically (as the duration of the busy period tends to infinity). Thus the behaviour of the number in the queue when it is non-empty is represented in [6], [7] by a probability density function $f(x,t)$ satisfying the diffusion equation

$$ (1) \qquad -\frac{\partial}{\partial t} f(x,t) - b \frac{\partial}{\partial x} f(x,t) + \frac{\alpha}{2} \frac{\partial^2}{\partial x^2} f(x,t) = 0 $$

where, $\{X(t), t \ge 0\}$ being the stochastic process approximating the number in queue,

$$ (2) \qquad f(x,t)dx = \Pr\{x \le X(t) < x+dx\} . $$

Since the approach was intended for heavy traffic conditions it was also assumed that the lower boundary at $x=0$ for the process $\{X(t), t \ge 0\}$ should act as a reflecting boundary. This last assumption implies that no probability mass can collect at $x=0$, so that the model becomes inadequate in predicting the probability of an empty queue. GAVER and
SHEDLER [8] and KOBAYASHI [9], who generalized this approach to queueing networks, have chosen the integration constants in the solution to (1) with reflecting boundaries so that the model correctly predicts the stationary probability of an empty queue. In [10] some results predicted by the diffusion approximation are compared with measurements obtained from simulation experiments and with known mathematical results illustrating the excellent degree of accuracy obtainable by this new approach.

In a previous paper [11] we proposed the use of FELLER'S elementary return process [12] instead of the diffusion process with reflecting boundaries to approximate a queueing system. This is a diffusion process with boundaries to which the process adheres for epochs which are exponentially distributed whenever the process attains a boundary; at the end of the epoch the process is reinitialized according to a fixed probability density function. This model proved to be self contained and capable of predicting a number of results already known from queueing theory (such as the ergodicity condition or the probability of an empty queue in stationary state) exactly. In approximating a queue, the epoch during which the process adheres to the boundary at $x=0$ (queue length is zero) is meant to represent an idle period, and this period is not exponentially distributed in general. This motivated the non-Markovian diffusion model studied in this paper.

In Sections 2 and 3 we develop a diffusion model which generalizes the model of FELLER by the introduction of COXIAN density functions for the holding times of the process at the boundaries, i.e. density functions whose Laplace-Stieltjes transform is a rational function [13]. We then apply this model to the approximation of the GI/G/1 queue and to a computer system model in Sections 4 and 5. In Section 6 we treat the case where the holding times have arbitrary continuous and differentiable density functions and prove an interesting and useful theorem: it is shown that the stationary probability distribution in this model depends only on the first moments of these densities.
## 2. A Non-Markovian Diffusion Process

The diffusion process we shall present in this section appears to be new, though it is a generalization of standard diffusion processes [12]. To simplify and motivate the description of our model we shall imagine that the stochastic process {$X(t), t \ge 0$}, which will be used later to approximate the number in queue, represents the position of a particle moving on the closed interval [0,M] of the real line. When the particle is in the open interval ]0,M[ its motion is described by a Wiener process with drift; b and $\alpha$, the instantaneous rate of change of the mean and variance of $X(t)$ in this region, are given by

$$ (3) \qquad b = \lim_{\Delta t \to 0} \frac{E\{X(t+\Delta t) - X(t)\}}{\Delta t} $$

$$ (4) \qquad \alpha = \lim_{\Delta t \to 0} \frac{E\{[X(t+\Delta t) - X(t)]^2\} - [E\{X(t+\Delta t) - X(t)\}]^2}{\Delta t} $$

For our present purposes it is not necessary that b and $\alpha$ be functions of x, t; this restriction can be relaxed however.

When the particle reaches the lower boundary of the interval [0,M] it remains there for a period of time h which is a random variable, at the end of which it "jumps" instantaneously back into the open interval ]0,M[ to a random point whose position is defined by the probability density function $f_l(x)$. The probability distribution function of h is

$$ (5) \qquad F_h(r) = \Pr\{h \le r\} $$

and is assumed to be a COX or almost general distribution function [13]. Since we shall be dealing with such distributions it may be useful to
recall here their definition and main properties; we shall do this by visualizing what happens to the particle which attains the boundary $x = 0$. We shall imagine that the particle hitting the lower boundary and staying there for a random time $h$ spends this time in traversing the network of Figure 1 which it enters at point A and leaves from B. This network is made up of a sequence of independent random delays each of which is exponentially distributed; the parameters of these exponential distribution functions are $\lambda_1, \dots, \lambda_n$, respectively. After traversing the i-th delay the particle enters the (i+1)-th delay with fixed probability $b_i$ or leaves the network of delays via B with probability $(1-b_i)$. Of course $b_n = 0$. The departure of the particle from the network of delays coincides with the end of $h$ the holding time at the lower boundary and the instantaneous jump of the particle back into the open interval $]0, M[$. Let us denote by $f_h(r)$ the probability density function of $h$. From the previous description we have that $f_h^*(s)$ the Laplace transform of $f_h(r)$ is given by

$$ (6) \qquad f_h^*(s) = \int_0^\infty e^{-sr} f_h(r) dr \\ = \sum_{i=1}^n (1-b_i) a_i \prod_{j=1}^i \frac{\lambda_j}{(s+\lambda_j)} $$

where

$$ (7) \qquad a_i = \begin{cases} 1 & \text{if } i = 1 \\ b_1 \cdots b_{i-1} & \text{if } i > 1 \end{cases} $$

In fact (6) is the partial fraction expansion of the Laplace transform of
a probability density function where the transform is a rational function of s if the $b_i$ and $\lambda_i$, $1 \le i \le n$, are allowed to be complex quantities. <fnref n="1" />⁾

Thus our *only* restriction on the random variable h is that the Laplace transform of its probability density function be a rational function.

Such density functions have been introduced by COX [14]. Their interest lies in the fact that they can be used to approximate a density function by matching its first K moments, where K can be arbitrarily large.

Also it is useful to note that there is a large body of literature on approximation by rational functions. Notice that the generality of (6) is limited by the fact that a positive probability that h is zero is not allowed; otherwise the expansion (6) would contain a constant term.

When the particle hits the upper boundary at $x = M$ it remains there for a random holding time H whose probability density function $f_H(r)$ is also almost general. The parameters in the network of delays corresponding to H, and similar to those for h shown in Figure 1, are $\mu_1, \dots, \mu_m$ and $B_1, \dots, B_m$ so that

$$ (8) \qquad \begin{aligned} f_H^*(s) &= \int_0^\infty e^{-sr} f_H(r) dr \\ &= \sum_{i=1}^m (1-B_i) A_i \prod_{j=1}^i \frac{\mu_j}{(s+\mu_j)} \end{aligned} $$

where

$$ (9) \qquad A_i = \begin{cases} 1 & \text{if } i = 1 \\ B_1 \cdots B_{i-1} & \text{if } i > 1. \end{cases} $$

<footnote n="1">¹⁾ As long as $f_h(r)$ is real.</footnote>
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 $f_2(x)$.

$f_1(x)$ and $f_2(x)$ may be taken to be functions of the instant at which the jumps occur. Notice that

$$ (10) \qquad E\{h\} = \sum_{i=1}^{n} (1-b_i) a_i \sum_{j=1}^{i} \frac{1}{\lambda_j} . $$

This may be rewritten as

$$ E\{h\} = \sum_{i=1}^{n} \frac{1}{\lambda_i} \sum_{j=i}^{n} (1-b_j) a_j $$

and it is easy to show that for $1 \le i \le n$,

$$ (11) \qquad a_i = \sum_{j=i}^{n} (1-b_j) a_j . $$

Therefore

$$ (12) \qquad E\{h\} = \sum_{i=1}^{n} \frac{a_i}{\lambda_i} . $$

Similarly

$$ (13) \qquad E\{H\} = \sum_{i=1}^{m} \frac{A_i}{\mu_i} . $$

Let us introduce the notation
$$ (14) \qquad \lambda = (E\{h\})^{-1}, \qquad \mu = (E\{H\})^{-1}. $$

If $\text{Var}(h)$ denotes the variance of the random variable $h$, we see that some simple algebra yields

$$ (15) \qquad \text{Var}(h) = \sum_{i=1}^{n} \frac{a_i}{\lambda_i^2}, \qquad \text{Var}(H) = \sum_{i=1}^{n} \frac{A_i}{\mu_i^2}. $$

That the process $\{X(t), t \ge 0\}$ defined in this section is non-Markovian is easy to see: once the particle is at any one of the boundaries the additional time it will remain there is not independent of the amount of time it has resided at the boundary up to the present instant.
Let $\{k(t), t \ge 0\}$ be the stochastic process such that for all $t \ge 0$,

$$ k(t) \in \{-n, -n+1, \dots, 0, 1, \dots, m\} $$

and

$$ (16) \quad k(t) = \begin{cases} -j, & 1 \le j \le n, \quad \text{if the particle is in the } j\text{-th stage of the holding time } h \\ 0, & \text{if the particle is in } ]0, M[ \\ i, & 1 \le i \le m, \quad \text{if the particle is in the } i\text{-th stage of the holding time } H. \end{cases} $$

Then the process $\{(X(t), k(t)), t \ge 0\}$ has the Markov property.
### 3. The Equations for the Model

Let $f = f(x,t)$ denote the probability density function of the stochastic process $\{X(t), t \ge 0\}$ in the open interval $]0,M[$ and let $A_{x,t}$ and $C_{x,t}$ be operators defined by

$$ (17) \qquad A_{x,t}f = - \frac{\partial}{\partial t} f - \frac{\partial}{\partial x} bf + \frac{1}{2} \frac{\partial^2}{\partial x^2} \alpha f $$

$$ (18) \qquad C_{x,t}f = -bf + \frac{1}{2} \frac{\partial}{\partial x} \alpha f . $$

Also let $P_i(t)$, $1 \le i \le 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 $Q_i(t)$, $1 \le i \le 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

$$ (19) \qquad A_{x,t}f + \sum_{i=1}^{n} \lambda_i (1-b_i) P_i(t) f_1(x) + \sum_{i=1}^{m} \mu_i (1-b_i) Q_i(t) f_2(x) = 0 $$

$$ (20) \qquad \frac{d}{dt} P_i(t) = \begin{cases} - \lambda_1 P_1(t) + C_{0,t} f, & \text{if } i = 1 \\ - \lambda_i P_i(t) + \lambda_{i-1} b_{i-1} P_{i-1}(t), & \text{if } 1 < i \le n. \end{cases} $$

$$ (21) \qquad \frac{d}{dt} Q_i(t) = \begin{cases} - \mu_1 Q_1(t) - C_M t f, & \text{if } i = 1 \\ - \mu_i Q_i(t) + \mu_{i-1} b_{i-1} Q_{i-1}(t), & \text{if } 1 < i \le m. \end{cases} $$

where

$$ C_{0,t}f = \lim_{x \to 0} \left[-bf + \frac{1}{2} \frac{\partial}{\partial x} \alpha f\right] $$

$$ C_M t f = \lim_{x \to M} \left[-bf + \frac{1}{2} \frac{\partial}{\partial x} \alpha f\right]. $$
Define $P(t)$ 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. We then have that

$$P(t) = \sum_{i=1}^{n} P_{i}(t), \quad Q(t) = \sum_{i=1}^{m} Q_{i}(t)$$

and

$$(22) \qquad \frac{d}{dt} P(t) = - \sum_{i=1}^{n} \lambda_{i} (1-b_{i}) P_{i}(t) + C_{0,tf}$$

$$(23) \qquad \frac{d}{dt} Q(t) = - \sum_{i=1}^{m} \mu_{i} (1-B_{i}) Q_{i}(t) - C_{M,tf}.$$

Equations (19), (20), (21) are simple to interpret. Suppose $\Omega$ is a subinterval of $]0, M[$. Then (19) can be deduced from

$$(24) \quad \frac{\partial}{\partial t} \int_{\Omega} f \, dx = \int_{\Omega} \left[ -\frac{\partial}{\partial x}bf + \frac{1}{2} \frac{\partial^2}{\partial x^2} \alpha f \right] dx + \sum_{i=1}^{n} \lambda_i (1-b_i) P_i(t) \int_{\Omega} f_1(x) dx \\ + \sum_{i=1}^{m} \mu_i (1-B_i) Q_i(t) \int_{\Omega} f_2(x) dx$$

which states that the rate of change of the probability mass in $\Omega$ is equal to the rate of flow of the probability mass via the boundaries of $\Omega$ (the first term on the right-hand-side of (24)) plus the rate of flow from $x=0$ and from $x=M$ (the second and third terms, respectively, on the right-hand-side). In order to deduce (20), notice that for $1 < i \le n$ we may write for any $t \ge 0$,
$$ (25) \quad P_i(t+\Delta t) = (1 - \lambda_i \Delta t) P_i(t) + \lambda_{i-1} \Delta t P_{i-1}(t) $$

since the time the particle spends in any one of the stages in Figure 1 is exponentially distributed; by collecting terms, dividing both sides by $\Delta t$ and taking $\Delta t \to 0$, (25) yields (20) for $1 < i \le n$ in the usual way. For obtaining (20), with $i=1$ a similar procedure is applied if one notices that $C_{0,t}f$ is the flow of probability mass out of $]0,M[$ from the lower boundary and, of course, into the first stage of the holding time at $x=0$. A similar interpretation can be given for (21); notice now that $-C_{M,t}f$ is the flow of probability mass away from $]0,M[$ via the upper boundary.

In addition to (19), (20), (21) appropriate boundary conditions for $f(x,t)$ must be specified and initial conditions (at $t=0$) must be given for the stochastic process. Since the boundaries at $x=0$ and $x=M$ behave as absorbing boundaries during their respective holding times we set $\lim_{x \to 0} f(x,t) = \lim_{M \to 0} f(x,t) = 0$ for all $t \ge 0$. Of course we set

$$ \int_0^M f \, dx + P(t) + Q(t) = 1 . $$

### 4. Application to the GI/G/1 Queue: Stationary Solution

In this section we propose an approximation to the number in a single server queue with general service time distribution of mean $1/\mu$ and variance $V_s$, independent of the interarrival times or of queue length, and with independent interarrival times having a general distribution function with mean $1/\lambda$ and variance $V_a$. The stochastic process
`{X(t), t ≥ 0}` approximating the number in queue at time `t` takes values on the non-negative real line `[0,∞[; it is of the type defined in the previous sections and its behavior is given by the equations

$$ (26) \qquad A_{x,t} f + \sum_{i=1}^{n} \lambda_i (1-b_i) P_i(t) \delta(x-1) = 0 $$

and

$$ (27) \quad \frac{d}{dt} P_i(t) = \begin{cases} - \lambda_1 P_1(t) + C_{0,t} f, & \text{if } i = 1 \\ - \lambda_i P_i(t) + \lambda_{i-1} b_{i-1} P_{i-1}(t), & \text{if } 1 < i \le n. \end{cases} $$

where $\delta(x)$ is the unit (probability) mass concentrated at $x=0$. In this model the case $X(t) = 0$ refers to the empty queue; an arrival at time $t$ to the empty queue corresponds to an instantaneous jump of $X(t)$ from 0 to 1, hence the value $f_1(x) = \delta(x-1)$. For a finite value of $t$ there can be no probability mass at infinity hence we only have a probability mass $P(t)$ at the origin. When the queue becomes empty, the general arrival process is being approximated by an almost general process in the sense of Section 2 of this paper with appropriately chosen parameters $\lambda_1, \dots, \lambda_n$ and $b_1, \dots, b_n$. The parameters $b$ and $\alpha$ in the operators $A_{x,t}$ and $C_{x,t}$ are chosen from the predictions of the central limit theorem as in [7,8,9]:

$$ (28) \qquad b = \lambda - \mu $$

$$ (29) \qquad \alpha = \lambda^3 v_a + \mu^3 v_s . $$

$\lambda$, the mean arrival rate to the queue, coincides with $(E\{h\})^{-1}$, where
h is the random variable representing the holding time at the lower boundary for the diffusion process. Since there is no upper boundary to the process being considered the definition of $\mu$ given in (14) is irrelevant to the discussion in this section. Let us denote by $K_a = \lambda^2 V_a$ and $K_s = \mu^2 V_s$ the squared coefficients of variation of the interarrival and service times, respectively.

Let us denote by $f(x)$ and $P_i$, $1 \le i \le n$, the stationary distribution for $\{X(t), t \ge 0\}$; it is obtained (when it exists) by setting $\frac{\partial f}{\partial t} = 0$, $\frac{dP_i(t)}{dt} = 0$, $1 \le i \le n$, in the system of equations (26), (27). The appropriate boundary condition is still $\lim_{x \to 0} f(x) = 0$. $P$ shall be the stationary total probability mass at the origin. The following result can be obtained.

**Proposition 1.** Let $\{X(t), t \ge 0\}$ be the stochastic process defined in this section approximating the number in queue for a GI/G/1 system. Its stationary probability distribution exists if $\rho = (\lambda/\mu) < 1$. Furthermore it is given by

$$ (30) \qquad f = \begin{cases} \rho[e^{-\gamma}-1]e^{\gamma x}, & \text{if } x \ge 1 \\ \rho[1-e^{\gamma x}], & \text{if } 0 \le x \le 1. \end{cases} $$

$$ (31) \qquad P = 1 - \rho $$

$$ (32) \qquad P_i = \frac{\mu a_i}{\lambda_i} \rho(1-\rho) = \frac{\lambda a_i}{\lambda_i} (1-\rho) $$

where
$$ (33) \qquad \gamma = - \frac{2(1-\rho)}{\rho K_a + K_s} . $$

The approximate expected queue length at stationary state is given by

$$ (34) \qquad \hat{L} = \int_{0}^{\infty} xf \, dx = \rho \left[ \frac{1}{2} + \frac{\rho K_a + K_s}{2(1-\rho)} \right] . $$

We shall not go through a proof of this result since it is very similar to that of Proposition 2 given in Section 5. It is interesting to notice that the diffusion approximation of this section is able to predict the stationary probability of an empty queue P exactly as one would expect to obtain from queueing theory; previous approaches to diffusion approximations [6,7,8,9] will not yield this result unless it is expressly introduced. Of course the fact that the total probability mass is unity

$$ P + \int_{0}^{\infty} f \, dx = 1 $$

has been used in obtaining our solution. Another interesting point is that $\hat{L}$ of (34) has a form similar to the POLLACZEK-KHINTCHINE formula for the M/G/1 queue which is

$$ L = \rho \left[ 1 + \frac{\rho(1+K_s)}{2(1-\rho)} \right] . $$

In fact, if we set $K_a = 1$ in (34) in order to approximate a Poisson arrival process we obtain that the error is
$$L - \hat{L} = \frac{\rho}{2}[1 - K_s] .$$

It is important to note that as $\rho \to 1$ the relative error $(L-\hat{L})/L$ vanishes, and that it vanishes as well when $K_s = 1$ (for the M/M/1 system). If we had discretized the probability density function $f(x)$ by defining the probability of $i$ customers in queue as

$$\pi_i = \int_{i-1}^{i} f dx, \quad \text{for } i \ge 1$$

and taken the average queue length as being

$$L' = \sum_{i=1}^{\infty} i \pi_i = \rho \left[ 1 + \frac{\rho K_s + K}{2(1-\rho)} \right]$$

we would have had

$$L - L' = - \frac{1}{2} \rho K_s .$$

Again the relative error vanishes as $\rho \to 1$.

What is most interesting is that the stationary probability of an empty queue and all the moments of the stationary queue length distribution predicted by Proposition 1 are independent of all moments but the first of the holding time distribution at the origin. 1)

It is essential to note that in our approximation method the random variable $h$ refers to the time interval between the last departure from the queue of a busy period to the first arrival of the next busy period. If the arrival process is Poisson it is natural to take $\lambda = (E\{h\})^{-1}$,

<footnote n="1">1) This will be proved in greater generality in Section 6.</footnote>
as we have done. 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 briefly the case where $(E\{h\})^{-1} = \lambda' \neq \lambda$. With this assumption the stationary solution of (26), (27) becomes

$$ (35) \qquad f = \begin{cases} R[e^{\lambda'}\gamma - 1]e^{\gamma x}, & x \geq 1 \\ R[1 - e^{\gamma x}], & 0 \leq x \leq 1 \end{cases} $$

$$ (36) \qquad P_i = \frac{\lambda'a_i}{\lambda_i}(1-R), \quad 1 \le i \le n $$

$$ (37) \qquad P = 1 - R $$

where

$$ (38) \qquad R = \frac{\lambda'}{\lambda' + \mu - \lambda} $$

The condition for existence of the stationary solution is still $\rho = \lambda/\mu < 1$ since it is derived from the condition $\gamma < 1$. We see however that the usual queueing theory result $P = 1 - \rho$ will only be obtained if we set $\lambda' = \lambda$.

## 5. Application to a Closed Two-Server System
### With General Service Time Distributions

A special case of the model presented in Section 3 will be proposed here as an approximation to a queueing system containing a finite number of customers and two servers.
The system whose behaviour we wish to approximate is shown in Figure 2. It consists of a central processing unit (CPU) and an input-output device (IOD); a finite fixed number M of programs are being executed in the system. We shall assume in this section that service times at the CPU are independent and identically distributed (i.i.d.) random variables with distribution function with mean $\mu^{-1}$ and variance $V_s$; they are independent of the service times at the IOD which are also i.i.d. random variables of mean $\lambda^{-1}$ and variance $V_a$. In general we do not exclude the possibility that $\lambda, V_a, \mu$ and $V_s$ be functions of M.

The approximate model is described by the following set of equations, where the operators $A_{x,t}$ and $C_{x,t}$ are given by (17) and (18).

$$ (39) \quad A_{x,t} f + \sum_{i=1}^{n} \lambda_i (1-b_i) P_i(t) \delta(x-1) + \sum_{i=1}^{m} \mu_i (1-B_i) Q_i(t) \delta(x-M+1) = 0 $$

$$ (40) \quad \frac{d}{dt} P_i(t) = \begin{cases} -\lambda_1 P_1(t) + C_{0,t} f, & \text{if } i=1 \\ -\lambda_i P_i(t) + \lambda_{i-1} b_{i-1} P_{i-1}(t), & \text{if } 1 < i \le n \end{cases} $$

$$ (41) \quad \frac{d}{dt} Q_i(t) = \begin{cases} -\mu_1 Q_1(t) - C_{M,t} f, & \text{if } i=1 \\ -\mu_i Q_i(t) + \mu_{i-1} B_{i-1} Q_{i-1}(t), & \text{if } 1 < i \le m \end{cases} $$

Here $f = f(x,t)$ is the probability density function for $\{X(t), t \ge 0\}$, the process approximating the number of programs in the CPU queue, when this number is in $]0,M[$. The epoch from the instant that the CPU queue
becomes empty to the arrival of the first customer is represented by the holding time h at the origin which is assumed to have an almost general distribution function (see Figure 1) of parameters $\lambda_1, \dots, \lambda_n$ and $b_1, \dots, b_n$. The corresponding parameters for H, the epoch from the instant that the IOD queue empties to the arrival of the first program to this queue, are $\mu_1, \dots, \mu_m$ and $B_1, \dots, B_m$. The boundary conditions are that $f(x,t)$ vanishes at $x=0$ and at $x=M$. As in the previous section, we take $b = \lambda - \mu$ and $\alpha = \lambda^3 v_a + \mu^3 v_s$ as the parameters in the operators $A_{x,t}$ and $C_{x,t}$.

**Proposition 2.** Let {$X(t), t \ge 0$} be the stochastic process approximating the number of programs in the CPU queue of the multiprogramming system described in this section, and whose probability distribution function satisfies (39), (40), (41) with the boundary conditions mentioned above. Its stationary probability density function exists and is given by

$$ (42) \qquad f = \begin{cases} K[1-e^{-\gamma x}], & 0 < x \le 1 \\ K[e^{-\gamma} - 1]e^{-\gamma x}, & 1 \le x \le M-1 \\ K[e^{\gamma(x-M)} - 1]e^{\gamma(M-1)}, & M-1 \le x < M \end{cases} $$

with P and Q, the probability masses at 0 and at M, respectively, at stationary state being

$$ (43) \qquad P = K \frac{(1 - \rho)}{\rho}, \quad Q = K (1 - \rho) e^{\gamma(M-1)} $$

where $\rho = \lambda / \mu$, and
$$K = \rho(1 - \rho^2 e^{\gamma(M-1)})^{-1}$$

This result can be either verified by substitution in (39), (40), (41) with the appropriate boundary conditions and setting partial derivatives with respect to time equal to zero or can be obtained by solving the differential equations directly. The verification of the result is left to the interested reader. It is interesting once again to notice that the stationary solution obtained depends only on the expected values of the holding times at 0 and M.

On Figure 3 is summarized the result of some simulation experiments conducted by M. BADEL and V.Y. SHUM concerning our model predictions for the model described on Figure 2. The quantity plotted is the absolute value of the error term relative to the quantity obtained by simulation, for the stationary probability (1-P) that the CPU is active. That is, if $η = (1-P)$ obtained from the diffusion model and $β$ is the corresponding quantity, then the quantity plotted is $|η - β|/β$. Two sets of simulation results, one with constant service time at the IOD and the other with exponentially distributed service time at the IOD are given. In each case we have also plotted the estimated confidence intervals for a 95% confidence level. The value of $ρ$ has been varied between 0.25 and 0.9 and M has been varied between 1 and 10; the relative error plotted for each value of $K_s$ is the maximum absolute relative error over all these values of $ρ$ and M for a given value of $K_s$. This error remains relatively low, and in any case is smaller than the width of the confidence interval.
## 6. The Diffusion Model with General Holding Times at the Boundaries

Consider the following equations :

$$ (44) \qquad \begin{aligned} \frac{\partial P(y,t)}{\partial t} &= \frac{\partial P(y,t)}{\partial y} + \frac{\partial P(y,t)}{\partial t} \\ &= -u(y) P(y,t) + [c_{o,t} f] \delta(y) \end{aligned} $$

$$ (45) \qquad \begin{aligned} \frac{\partial Q(z,t)}{\partial t} &= \frac{\partial Q(z,t)}{\partial z} + \frac{\partial Q(z,t)}{\partial t} \\ &= -v(z) Q(z,t) - [c_{M,t} f] \delta(z-M) \end{aligned} $$

$$ (46) \qquad A_{x,t} f + \left( \int_{0}^{\infty} u(y) P(y,t) dy \right) f_1(x) + \left( \int_{0}^{\infty} v(z) Q(z,t) dz \right) f_2(x) = 0 $$

where $f(c,t) = f(M,t) = 0$, and

$$ \int_{0}^{M} f(x,t) dx + \int_{0}^{\infty} P(y,t) dy + \int_{0}^{\infty} Q(z,t) dz = 1 $$

$$ (47) \qquad u(y) = \frac{g(y)}{1 - G(y)}, \quad v(z) = \frac{h(z)}{1 - H(z)} $$

$G(y)$, $H(z)$ are the probability distribution functions of total time spent by the process $\{X_t, t \ge 0\}$ on the lower and upper boundary, respectively, and $g(y)$, $h(z)$ are the corresponding densities (which are assumed to exist). The above system of equations generalizes (19), (20), (21) to the case where the holding time probability distribution functions $G(y)$, $H(z)$ are continuous but not necessarily Coxian. $f_1(x)$ and $f_2(x)$ are the density functions (defined on $[0, M]$) for the value taken by $\{X_t, t \ge 0\}$ at the time just after a jump from the lower and upper boundaries, respectively. The probabilities $P(y,t)$, $Q(z,t)$ are:

$$ P(y,t) = \text{Pr} [X_{\tau} = 0, \tau \in ]t - y, t] \text{ and } X_{t - y} > 0] $$

$$ Q(z,t) = \text{Pr} [X_{\tau} = M, \tau \in ]t - z, t] \text{ and } X_{t - z} < M] $$
The following result states that the stationary solution to (44), (45), (46) depends only on the expected values of the holding times at the upper and lower boundaries. Let

$$ 
\begin{aligned}
\lambda^{-1} &= \int_{0}^{\infty} t g(t) \, dt, \quad \mu^{-1} = \int_{0}^{\infty} t h(t) \, dt \\
P(t) &= \int_{0}^{\infty} dy P(y, t), \quad Q(t) = \int_{0}^{\infty} dz Q(z, t)
\end{aligned}
$$

$P(t)$, $Q(t)$ are the probabilities that $X_t = 0$, $X_t = M$, respectively.

**Theorem** : Consider the stationary solution to (44), (45), (46) for $P(t)$, $Q(t)$, $f(x,t)$, obtained by setting

$$ \frac{\partial P(y, t)}{\partial t} = 0, \quad \frac{\partial Q(z, t)}{\partial t} = 0, \quad \frac{\partial f}{\partial t} = 0 $$

It depends on the expected values of the holding times at $x=0$, $x=M$, and is independent of all higher moments of the holding time distributions. Thus the stationary solution is identical to the corresponding solution when the holding times are exponentially distributed with means $\lambda^{-1}$, $\mu^{-1}$.

**Proof.** : The stationary solution to (44) is $P(y) = \lim_{t \to \infty} P(y,t) = a(y)\exp(-\int_0^y u(s)ds)$ yielding $P(y) = a(y)(1-G(y))$, where $a(y)$ is obtained using the boundary condition $f(0,t)=0$ in (44) and carrying out the usual substitution, so that

$$ a(y) = \left[ \frac{\partial}{\partial x} (\alpha f) \right]_x = 0 $$

so that

$$ P = \int_{0}^{\infty} P(y) \, dy = \frac{1}{2} \left[ \frac{\partial}{\partial x} (\alpha f) \right]_{x=0} $$

similarly

$$ Q = \frac{1}{2} \left[ - \frac{\partial}{\partial x} (\alpha f) \right]_x = M $$

also we obtain

$$ \int_{0}^{\infty} u(y) P(y) \, dy = \frac{1}{2} \left[ \frac{\partial}{\partial x} (\alpha f) \right]_{x=0} = \lambda P $$

$$ \int_{0}^{\infty} v(z) Q(z) \, dz = \mu Q $$

so that (46) becomes at stationary state

$$ (48) \quad - \frac{\partial(bf)}{\partial x} + \frac{\partial^2(a f)}{\partial x^2} + \lambda P f_1(x) + \mu Q f_2(x) = 0 $$
which is independent of all moments but the first of $G(y)$, $H(z)$. This completes the proof.

This result is particularly interesting since it states that if we are only interested in the equilibrium solution to the non-Markovian diffusion model it suffices to solve the corresponding Markovian diffusion model.
<figure><img src="image_3.png" /><figcaption>Figure 1</figcaption></figure>

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

1. Jackson, J.R. "Jobshop-like queueing systems," Management Science 10, 131-142, 1963.

2. Moore, C.G. "Queueing network models of time-sharing systems," Ph.D. Thesis, Computer Science Dept., University of Michigan, Ann Arbor, 1970.

3. Muntz, R.R. "Poisson departure processes and queueing networks," IBM Research Report RC 4145, December 1972.

4. Chandy, K.M. "The analysis and solutions for general queueing networks," Proc. Sixth Annual Princeton Conf. on Information Sciences and Systems, Princeton University, Princeton, 1972.

5. Baskett, F., Chandy, K.M., Muntz, R.R., Palacios, J. "Open, closed, and mixed networks of queues with different classes of customers," JACM 22, 2, April 1975.

6. Gaver, D.P. "Diffusion approximations and models for certain congestion problems," J. Appl. Prob. 5, 607-623, 1968.

7. Newell, G.F. "Applications of queueing theory," Chapman and Hall, London, 1971.

8. Gaver, D.P., Shedler, G.S. "Processor utilization in multiprogramming systems via diffusion approximations," Operations Research 21, 569-576, 1963.

9. Kobayashi, H. "Application of the diffusion approximation to queueing networks: Part I- Equilibrium queue distributions," JACM 21, 2, 316-328, April 1974.

10. Reiser, H., Kobayashi, H. "Accuracy of a diffusion approximation for some queueing networks," IBM J. Res. Dev. 18, 2, 110-124, March 1974.

11. Gelenbe, E. "On approximate computer system models," JACM 22, 2, 261-269, April 1975.

12. Feller, W. "Diffusion processes in one dimension," Trans. Am. Math. Soc., 77, 1-31, 1954.

13. Cox, D.R. "A use of complex probabilities in the theory of stochastic processes," Proc. Cambridge Philosophical Society, 51, 2, 313-319, April 1955.

14. Gelenbe, E., Muntz, R.R. "Probabilistic models of computer systems: Part I (Exact results)," to appear in Acta Informatica.

15. Gelenbe, E., Pujolle, G. "Probabilistic models of computer systems: Part II (Approximating the behaviour of a single queue in a network)," to appear in Acta Informatica.

16. Gelenbe, E. "A diffusion model for drum input-output operations with general interarrival distribution," in preparation.