An approximate numerical solution for multiclass preemptive priority queues with general service time distributions
Authors: Snyder Patricia M. , Stewart William J. 
Published: January 2007

In this paper an approximate numerical solution for a multiclass preemptive priority single server queue is developed. The arrival process of each class follows a Poisson distribution. The service time distribution must have a rational Laplace transform, but is otherwise arbitrary and may be different for different classes. The work reported here was motivated by a desire to compute the equilibrium probability distribution of networks containing preemptive priority servers. Such networks are frequently encountered when modeling computer systems, medical care delivery systems and communication networks. We wish to use an iterative technique which constructs a series of two station networks consisting of one station from the original network and one “complementary” station whose behavior with respect to the original station mimics that of the rest of the network. At each iteration, it is necessary to compute the equilibrium probability distribution of one or more preemptive priority queues.
Although such queues have been studied for some time, the resulting solutions have most often been developed utilizing transforms or probability generating functions, e.g. Jaiswal [1968]; in many cases of interest, inversion has not been attempted. Miller [1981] presented explicit solutions for two class priority queues but Miller's work, which is based on that of Neuts, is limited to exponential service times and two classes. The approach presented here is applicable to many classes and to more general service time distributions than have previously been considered.

The algorithm utilizes a bootstrap approach, a concept borrowed from dynamic programming. The solution for class 1 is trivial. Once we have solved the system with k different classes, we have available all the necessary information to solve the system with k+l classes. We shall assume that each class has a distinct service time distribution G
k
, with mean g
k
and variance s
2
k
. Let class k have preemptive priority over class l if k < l.

Successive steps of the algorithm are based upon the machine breakdown and repair model, previously used by Keilson [1962]; White and Christie [1958]; Gaver [1962]; and Avi-Itzhak and Noar [1963] to model preemptive priority queues. When we are solving for the equilibrium probability distribution of class k jobs, k > 1, we consider a model with one machine whose service time distribution is G
k
. The breakdown rate of the machine is the sum of the arrival rates of all higher priority jobs. The downtime or repairtime of the machine has mean γ
k-l
and variance σ
2
k-l
; these parameters are the mean and variance of the busy period in a preemptive priority system with k-l classes.

The first step in the solution of all the machine breakdown and repair models considered herein is to construct the infinitesimal generator @@@@ which specifies the transitions into and out of each state for the given model. Given @@@@, it is a simple matter to write the global balance equations, i.e., the Chapman-Kolmogorov equations which relate flow into and out of each state of the model at equilibrium. These global balance equations are second order difference equations, i.e. if @@@@
k
(n) is the probability of finding n class k jobs in the system at equilibrium, @@@@
k
(n) is a function of @@@@
k
(n-1) and @@@@
k
(n+1). The next step is to construct another set of balance equations, which we call the aggregate balance equations, and to use the results in Snyder and Stewart [1985] to reduce each of the second order difference equations (the global balance equations) to first order difference equations. As long as the generator @@@@ is irreducible and positive recurrent, it is now possible to define a coefficient matrix @@@@
k
which relates the probability vectors @@@@
k
(n) and @@@@
k
(n-1); specifically, we shall construct real matrices @@@@
k
whose elements are given explicitly as a function of the model parameters and for which @@@@
k
(n) = @@@@
k
(n-1) @@@@
k
.