0

Zuk and Kirszenblat

Non-Preemptive Priority Queue

Explicit Results for the Distributions of Queue Lengths for a Non-Preemptive
Two-Level Priority Queue

Josef Zuk Defence Science and Technology Group, Melbourne, Australia,
<EMAIL_ADDRESS>David Kirszenblat Defence Science and Technology
Group, Melbourne, Australia<EMAIL_ADDRESS>
Explicit results are derived using simple and exact methods for the joint and
marginal queue-length distributions for the M/M/$c$ queue with two non-
preemptive priority levels. Equal service rates are assumed. Two approaches
are considered. One is based on numerically robust quadratic recurrence
relations. The other is based on a complex contour-integral representation
that yields exact closed-form analytical expressions, not hitherto available
in the literature, that can also be evaluated numerically with very high
accuracy.

queueing theory; non-preemptive priority; queue length distribution Primary:
90B22; secondary: 60K25, 60J74 Primary: Queues: Priority; secondary: Queues:
Markovian

Date created: July 03, 2023. Last update: September 14, 2023.

##  1 Introduction

This work is concerned with the development of practical algorithms for the
computation of joint and marginal distributions of queue lengths for the
M/M/$c$ queue with a non-preemptive priority discipline. Applications of this
model are found in telecommunications [5], health care [10, 25], radar [21,
22], air traffic control [23] and numerous other areas.

The non-preemptive priority queue discipline is as stated by Dressin and Reich
[7]: Once a client’s service has begun, it is permitted to proceed to
completion. If a server becomes empty, and there is at least one client
waiting, then a client of the highest priority present in the queue is
admitted to the server. Clients of equal priority are served on a first-come,
first-served basis. This is also known as the ‘head of the line’ discipline.

Thus, let us consider a non-preemptive queue with K priority levels, each with
a distinct Poisson arrival rate $\lambda_{k}$, $k=1,2,\ldots,K$ and
corresponding level traffic intensity111Consistent with [8, 11], the notation
$\rho_{k}$ reserved for $\rho_{k}\equiv\lambda_{k}/\mu$, so that
$r_{k}=\rho_{k}/N$. $r_{k}=\lambda_{k}/(N\mu)$, leading to a total traffic
intensity for the aggregation of all arrivals of $r=\sum_{k=0}^{K}r_{k}$. We
adopt the usual convention that smaller priority-level indices $k$ represent
higher priorities. Thus, $r_{1}$ denotes the traffic intensity associated with
the highest priority level. For simplicity, we have assumed a common
exponential service rate $\mu$ among all priority levels. The number of
servers is denoted by $c=N$.

In this work, attention is confined to the two-level problem $K=2$. Analysis
of this case is amenable to a number of analytical techniques that do not
extend easily, or at all, to the general multi-level priority problem. Also,
the two-level problem has a distinguished status. All marginal distributions
for the multi-level problem can be inferred from the low-priority marginal
pertaining to just two priority levels [6]. If we let $r_{\text{hi}}$ and
$r_{\text{lo}}$ denote the level traffic intensities for the high and low
priority arrivals, respectively, for the two-level problem, then the wait-
conditional222See below in Section 2. marginal distribution of the queue
length for priority level $k=1,2,\ldots,K$ in the multi-level problem is
obtained by making the identifications

$r_{\text{lo}}=r_{k}\;,\quad r_{\text{hi}}=\sum_{\ell=1}^{k-1}r_{\ell}\;,$ (1)

so that the total traffic intensity in the effective (wait-conditional) two-
level problem becomes $r=r_{\text{sum}}$, with

$r_{\text{sum}}=r_{\text{lo}}+r_{\text{hi}}=\sum_{\ell=1}^{k}r_{\ell}\;.$ (2)

For the actual two-level problem, we have the identifications
$r_{\text{hi}}\equiv r_{1}$, $r_{\text{lo}}\equiv r_{2}$, and we shall use
both sets of notation interchangeably. It is also convenient to introduce the
parameter $\nu$ that represents the fraction of all arrivals that are of high
priority (which we abbreviate as ‘hifrac’). Thus $r_{\text{hi}}=\nu r$,
$r_{\text{lo}}=(1-\nu)r$, $0\leq\nu\leq 1$.

Previous work on the non-preemptive priority queue has focused, almost
entirely, on calculating moments and the waiting-time distributions per
priority level. In early work, Cobham [3, 4], followed by Holley [9], were
first to consider the mean waiting times and queue lengths. Waiting-time means
and second moments for general service-time distributions were subsequently
given by Kesten and Runnenberg [15]. Gail et al. [8] studied the non-
preemptive M/M/$c$ system for two priority levels with different exponential
service rates. While they developed a matrix algorithm for determining various
characteristics of a generating function for this problem, explicit results
were also confined to the mean waiting times and queue lengths.

For the waiting-time problem, Davis [6] improved on previous work by Dressin
and Reich [7] to derive an explicit integral expression for the waiting-time
distribution for the non-preemptive priority queue. He analysed the two-level
problem, as the waiting time distribution for the multi-level problem can be
inferred from the two-level case. He did not study the queue-length marginals,
and they cannot be directly inferred from the waiting-time distributions by
appealing to the distributional form of Little’s law [1, 13] as the no-
overtaking assumption is violated. Kella and Yechiali [14] covered the same
ground as Davis for the probability waiting function (PGF) of the waitingtime,
but using a different methodology. The moment generating function (MGF) of the
waiting time and associated moments have also been considered in [19]. More
recently, Wagner [26] has studied the waiting-time MGF for a finite-capacity,
multi-server version of the same problem as Davis.

For the queue-length distributions, Miller [17, 18] uses a matrix-geometric
method for the two-level problem that results in a complex algorithm involving
multiple levels of recursion. Little is presented about the numerical
stability of this approach, and it is known to deteriorate for traffic
intensities close to unity. Kao and Narayanan [11] and Kao and Wilson [12]
also employ matrix-geometric methods for the two-level problem which, as they
point out, unavoidably require finite-state truncation. The aforementioned
papers deal with unequal service rates. The matrix-geometric method [20]
applied to queueing models has the singular disadvantage that it necessitates
truncation of the problem to prescribed finite maximum values of queue
lengths. While powerful, it is complex and not elegant. Thus, its use should
best be avoided whenever simpler alternatives are available, and this is
manifestly the case for the present problem, as will become clear.

In earlier work, Marks [16] studied the two-level problem with common service
rate and derived a highly complex system of linear partial difference
equations that must be solved recursively. The required manipulations are
cumbersome and no insight into the analytic structure of the problem is
gained. However, it is most likely the first paper where actual queue-length
probabilities, rather than the PGF, were computed. No light is shed on the
numerical stability of the method.

A different approach, based on a partial PGF, is due to Cohen [5], who studied
the two-level problem with equal service rates; and it is this approach that
we pursue in the discussion that follows. We take up the programme where Cohen
[5] left off, in devising simple and practical schemes for extracting actual
probabilities from the PGF. There is the additional benefit that this approach
can be extended to the general multi-level problem. Shortle et al. [24] have
remarked that ‘the determination of stationary probabilities in a non-
preemptive Markovian system is an exceedingly difficult matter, well near
impossible when the number of priorities exceeds two’. In a separate
forthcoming paper, we shall demonstrate otherwise.

The present work focuses on explicit results that are useful for practical
applications. While we do not purport to have made general theoretical
advances in priority queues, the work does serve to fill a large gap in the
literature by establishing basic results for a paradigmatic model that one
would expect to have been uncovered decades ago. We believe that it also has
pedagogical value. For the two-level non-preemptive priority queue, Shortle et
al. [24], in the most recent edition of their textbook, set up the stationary
balance equations but remark that ‘obtaining a reasonable solution to these
stationary equations is very difficult, … The most we can do comfortably is
obtain expected values via two-dimensional generating functions’. The
simplicity of the methods described herein might render a more detailed
treatment of the subject suitable for elementary texts.

##  2 Non-Preemptive Priority Queue

The no-wait probability $P_{\text{NW}}$ is the probability that a new arrival
will find at least one server idle. It is clearly independent of the queue
discipline, and is given by [6]

$\frac{1}{1-P_{\text{NW}}}=1+(1-r)\frac{N!}{(Nr)^{N}}{\cdot}\sum_{k=0}^{N-1}\frac{(Nr)^{k}}{k!}\;.$
(3)

Let $P(n,m)$ denote the steady-state probability that there are $n$ low-
priority clients in the queue (rather than in the system) and $m$ high-
priority clients in the queue. We have the decomposition [cf. 6]

$P(n,m)=P_{\text{NW}}{\cdot}\delta_{n0}\delta_{m0}+(1-P_{\text{NW}}){\cdot}f(n,m)\;,$
(4)

where $f(n,m)$ represents the wait-conditional joint PMF, i.e. the probability
that there are $n$ low-priority clients and $m$ high-priority clients in the
queue, given that all servers are busy. The wait-conditional distribution does
not explicitly depend on the number of servers $N$. There is only an indirect
dependence on $N$ through the total traffic intensity $r$.

Our starting point is the paper of Cohen [5], which introduced a partial PGF
for the problem that summed only over the low-priority argument:

$G_{m}(p)\equiv\sum_{n=0}^{\infty}p^{n}P(n,m)\;.$ (5)

This turns out to be a very convenient strategy, especially given the fact
that the wait-conditional high-priority marginal is a simple geometric
distribution. Only the low-priority marginal is non-trivial. We introduce a
wait-conditional version $g_{m}(p)$ of this PGF such that

$g_{m}(p)\equiv\sum_{n=0}^{\infty}p^{n}f(n,m)\;.$ (6)

It follows that

$G_{m}(p)=P_{\text{NW}}{\cdot}\delta_{m0}+(1-P_{\text{NW}}){\cdot}g_{m}(p)\;.$
(7)

Cohen’s result [5] for the wait-conditional PGF for the two-level non-
preemptive priority queue with equal service rates is333The quantities
$\lambda_{1,2}(p)$ should not be confused with the arrival rates introduced
earlier. We are adhering to Cohen’s original, but less than ideal, notation.

$g_{m}(p)=\frac{(1-r)(1-p)}{1-p\lambda_{2}(p)}{\cdot}\lambda_{1}^{m}(p)\;,$
(8)

where $\lambda_{1,2}(p)$ are defined as follows: Let us introduce
$\lambda(p)=\lambda_{\pm}(p)$ as the two solutions of the quadratic equation

$\lambda^{2}-(1+r-r_{2}p)\lambda+r_{1}=0\;,$ (9)

such that

$\lambda_{\pm}(p)=[b(p)\pm\sqrt{b^{2}(p)-4r_{1}}]/2\;,\quad b(p)\equiv
1+r-r_{2}p\;.$ (10)

Then in (8), we have $\lambda_{1}(p)=\lambda_{-}(p)$,
$\lambda_{2}(p)=\lambda_{+}(p)$, and it is useful to note that

$\lambda_{+}(p)+\lambda_{-}(p)=b(p)\;,\quad\lambda_{+}(p)\cdot\lambda_{-}(p)=r_{1}\;.$
(11)

Another way to express the PGF for the wait-conditional distribution is

$g_{m}(p)=g_{\text{lo}}(p){\cdot}[1-\lambda_{1}(p)]\lambda_{1}^{m}(p)\;,\quad
g_{\text{lo}}(p)=\frac{1-r}{\lambda_{2}(p)-r}\;.$ (12)

That $g_{\text{lo}}(p)$ represents the wait conditional PGF for the low-
priority marginal is clear from observing that

$\sum_{m=0}^{\infty}g_{m}(p)=g_{\text{lo}}(p)\;.$ (13)

On the other hand, it follows directly from (12) that the wait-conditional PMF
for the high priority marginal is given by

$f_{\text{hi}}(m)=g_{m}(1)=(1-r_{1})r_{1}^{m}\;.$ (14)

Consequently, the only marginal distribution of interest in the present study
is that for the low-priority level.

By construction, the wait-conditional joint PMF $f(n,m)$ is recovered from the
PGF $g_{m}(p)$ according to

$f(n,m)=\frac{1}{n!}{\cdot}\left.\frac{d^{n}}{dp^{n}}g_{m}(p)\right|_{p=0}\;.$
(15)

The multiple derivative is prohibitively cumbersome to directly perform
analytically. Thus, we proceed to present two alternative strategies that
render the problem tractable.

##  3 Quadratic Recurrence

The first method constructs a recurrence relation based on the fact that the
functions $\lambda_{\pm}(p)$ solve a quadratic equation. We begin by
considering the low-priority marginal, whose PGF can be expressed as

$g_{\text{lo}}(p)=\frac{1-r}{\lambda_{2}(p)-r}\;.$ (16)

Since $\lambda_{2}(p)$ satisfies a quadratic equation, then so does
$g_{\text{lo}}(p)$. Let us set

$u\equiv r_{2}p\;,\quad
g(u)\equiv\frac{1}{\lambda_{2}-r}=\sum_{k=0}^{\infty}\frac{g_{k}}{k!}u^{k}\quad\Rightarrow\quad
g_{k}=\left.\frac{d^{k}g(u)}{du^{k}}\right|_{u=0}\;.$ (17)

Then we obtain

$(ru-r_{2})g^{2}+(u-1+r)g+1=0\;.$ (18)

We now differentiate this equation $n$ times with respect to $u$, and use the
identities

$\displaystyle\begin{aligned}
\frac{1}{n!}{\cdot}\left.\frac{d^{n}}{du^{n}}(ug)\right|_{u=0}&=\frac{g_{n-1}}{(n-1)!}\;,\\\
\frac{1}{n!}{\cdot}\left.\frac{d^{n}}{du^{n}}(g^{2})\right|_{u=0}&=\sum_{k=0}^{n}\frac{g_{k}}{k!}{\cdot}\frac{g_{n-k}}{(n-k)!}\;,\\\
\frac{1}{n!}{\cdot}\left.\frac{d^{n}}{du^{n}}(ug^{2})\right|_{u=0}&=\sum_{k=0}^{n-1}\frac{g_{k}}{k!}{\cdot}\frac{g_{n-k-1}}{(n-k-1)!}\;.\end{aligned}$
(19)

For the quantities $f_{k}\equiv g_{k}/k!$, this leads to the non-linear
recurrence relations

$f_{n}=\frac{1+rf_{0}}{1-r+2r_{2}f_{0}}{\cdot}f_{n-1}+\frac{1}{1-r+2r_{2}f_{0}}\sum_{k=1}^{n-1}f_{k}{\cdot}(rf_{n-k-1}-r_{2}f_{n-k})\;,$
(20)

for $n=1,2,\ldots$, with

$f_{0}=\frac{1}{2r_{2}}\left[\sqrt{(1-r)^{2}+4r_{2}}-(1-r)\right]=\frac{2}{1-r+\sqrt{(1-r)^{2}+4r_{2}}}>0\;.$
(21)

The expression for $f_{0}$ follows from
$f_{0}^{-1}=g_{0}^{-1}=\lambda_{2}(0)-r$. We observe that
$f_{\text{lo}}(n)=(1-r)r_{2}^{n}f_{n}$.

Efficient vectorized implementations in Matlab are possible. Practical
implementation proceeds as follows: Let us introduce an arbitrary scale factor
$\Lambda$, define

$D\equiv 1-r+2r_{2}f_{0}=\sqrt{(1-r)^{2}+4r_{2}}\;,\quad c_{1}\equiv
r_{2}/\Lambda\;,\quad c_{2}\equiv\Lambda/D\;,$ (22)

and scale according to
$\tilde{f}_{n}\equiv\Lambda^{n}f_{n}=(r_{2}/c_{1})^{n}f_{n}$. Then we solve
the recurrence

$\tilde{f}_{n}=c_{2}{\cdot}\biggl{(}\tilde{f}_{n-1}+\sum_{k=0}^{n-1}\tilde{f}_{k}\Delta_{(n-1)-k}\biggr{)}$
(23)

and recover the marginal as $f_{\text{lo}}(n)=(1-r)c_{1}^{n}\tilde{f}_{n}$. At
each step, we set

$\Delta_{k}\equiv r\tilde{f}_{k}-c_{1}\tilde{f}_{k+1}\;,$ (24)

for $k=1,2,\ldots,n-1$, subject to the initialization $\tilde{f}_{n}\leftarrow
0$ within the scope of evaluating $\Delta_{n-1}$. We find that good numerical
performance is achieved with $\Lambda=r_{2}$, so that $c_{1}=1$.

Analogous treatment of the joint PMF is only marginally more complex. Based on
the quadratic

$\lambda_{\pm}^{2}+(u-1-r)\lambda_{\pm}+r_{1}=0\;,$ (25)

we solve for the Taylor-series coefficients $\lambda^{(k)}_{\pm}$ in

$\lambda_{\pm}=\sum_{k=0}^{\infty}\lambda^{(k)}_{\pm}p^{k}=\sum_{k=0}^{\infty}\Lambda^{-k}f^{\pm}_{k}u^{k}\;,$
(26)

for some arbitrary scale factor $\Lambda$, using the non-linear recurrence

$f^{\pm}_{n}=\mp\frac{1}{\sqrt{(1-r)^{2}+4r_{2}}}\left(\Lambda
f^{\pm}_{n-1}+\sum_{k=1}^{n-1}f^{\pm}_{k}{\cdot}f^{\pm}_{n-k}\right)\;,$ (27)

$n=1,2,\ldots$, where

$f^{\pm}_{0}=\tfrac{1}{2}\left(1+r\pm\sqrt{(1-r)^{2}+4r_{2}}\right)\;.$ (28)

The $\lambda$-coefficients are recovered according to
$\lambda^{(k)}_{\pm}=(r_{2}/\Lambda)^{k}f^{\pm}_{k}$. As with the marginal,
the choice $\Lambda=r_{2}$ results in good numerical performance. All that
remains to be done is to use the standard recursion for multiplication of
power series as dictated by (12). The simplest way to proceed is via repeated
convolutions:

$\displaystyle\begin{aligned} \phi_{0}&=(1-r){\cdot}\mathop{\rm
conv}\nolimits\left(\frac{1}{\lambda_{2}-r},1-\lambda_{1}\right)\;,\\\
\phi_{k}&=\mathop{\rm conv}\nolimits(\phi_{k-1},\lambda_{1})\;,\end{aligned}$
(29)

for $k=1,2,\ldots,m$. Then $f(n,m)=\phi_{m}(n)$. The conv function is defined
like the Matlab function of the same name: Suppose that $C(u)=A(u)B(u)$, with

$A(u)=\sum_{n=0}^{n_{1}}a(n)u^{n}\;,\quad
B(u)=\sum_{n=0}^{n_{2}}b(n)u^{n}\;,\quad
C(u)=\sum_{n=0}^{n_{1}+n_{2}}c(n)u^{n}\;.$ (30)

Then $c=\mathop{\rm conv}\nolimits(a,b)$, where

$c(n)=\mathop{\rm conv}\nolimits(a,b)(n)\equiv\sum_{k=0}^{n}a(k)b(n-k)\;,$
(31)

for $n=0,1,2,\dots,n_{1}+n_{2}$.

The foregoing recurrence relations constitute a significant improvement over
the strategy implemented in [2], and are vastly simpler than those arising
from the matrix-geometric method as considered in [11, 12, 17, 18]. While the
quadratic recurrence method exhibits excellent numerical behaviour, it gives
little insight into the analytical structure of the distributions. This
deficiency is addressed in the next section.

##  4 Complex Contour Integral

x The $z$-contour that results from taking the $p$-contour to be the unit
circle centred on the origin, plotted for the case of total traffic intensity
$r=0.95$ and fraction of high-priority arrivals $\nu=0.75$. Also displayed are
the locations of R-integral poles $z_{0},z_{1},z_{2}$.

Another strategy in dealing with (15) is to represent it in terms of a complex
contour integral in accordance with Cauchy’s integral theorem. This yields

$f(n,m)=(1-r)\oint_{\mathcal{C}}\frac{dp}{2\pi
i}\,\frac{(1-p)\lambda_{1}^{m}}{p^{n+1}(1-p\lambda_{2})}\;,$ (32)

where $\mathcal{C}$ is an anti-clockwise circle centred about the origin with
radius less than $1/r$. It follows directly that the low-priority marginal
PMF, defined by

$f_{\text{lo}}(n)\equiv\sum_{m=0}^{\infty}f(n,m)\;,$ (33)

is represented as a complex contour integral by

$f_{\text{lo}}(n)=(1-r)\oint_{\mathcal{C}}\frac{dp}{2\pi
i}\,\frac{1-p}{p^{n+1}}{\cdot}\frac{1}{(1-p\lambda_{2})(1-\lambda_{1})}\;.$
(34)

The conventional approach in dealing with such contour integrals, mirroring
the approach adopted previously for the waiting-time distribution [6], would
be to deform the contour by expanding it to the circle at infinity while
avoiding a cut of finite extent on the real axis that is generated by the
square-root component of $\lambda_{\pm}(p)$, and a possible simple pole that
also lies on the real axis. The circle at infinity yields a vanishing
contribution, which leaves a (potential) pole term and a real-valued integral
along the cut. We shall explore this approach separately in a forthcoming
paper, where we shall show that it leads to integral expressions that are
amenable to efficient quadrature algorithms, and can also be evaluated
analytically in terms of a generalized form of the associated Legendre
functions. In the present work, we pursue a different method based on a change
of integration variable.

Let $\lambda=z_{1},z_{2}$ be the roots of the polynomial equation
$\lambda^{2}-(1+r)\lambda+r_{1}=0$, so that we have $z_{1}+z_{2}=1+r$,
$z_{1}z_{2}=r_{1}=\nu r$. Then, the inversion of $z=\lambda_{\pm}(p)$ yields

$p=-(z-z_{1})(z-z_{2})/(r_{2}z)\;.$ (35)

Thus,

$dp=-\frac{1}{r_{2}}\left(1-\frac{z_{1}z_{2}}{z^{2}}\right)dz\;,$ (36)

in which case

$\frac{dp}{p^{n+1}}=(-r_{2})^{n}\frac{(z^{2}-z_{1}z_{2})z^{n-1}}{\left[(z-z_{1})(z-z_{2})\right]^{n+1}}{\cdot}dz\;.$
(37)

We make the change of integration variable $p\mapsto z:z=\lambda_{1}(p)$, in
which case $\lambda_{2}(p)=r/z$, and we make the identifications

$z_{0}=r_{1}/r\;,\quad z_{1}=\lambda_{-}(p=0)\;,\quad
z_{2}=\lambda_{+}(p=0)\;,$ (38)

or, equivalently,

$z_{0}=\nu\;,\quad z_{1}=\tfrac{1}{2}\left[1+r-\sqrt{(1+r)^{2}-4\nu
r}\right]\;,\quad z_{2}=\tfrac{1}{2}\left[1+r+\sqrt{(1+r)^{2}-4\nu
r}\right]\;.$ (39)

Then, we obtain

$\frac{1-p}{1-p\lambda_{2}}=\frac{z}{r}{\cdot}\frac{z-1}{z-z_{0}}\;.$ (40)

It follows that the joint PMF is given by

$f(n,m)=\frac{(1-r)(-r_{2})^{n}}{r}\oint_{\mathcal{C}^{\prime}}\frac{dz}{2\pi
i}\,\frac{z^{m+n}}{z-z_{0}}{\cdot}\frac{(z-1)(z^{2}-z_{1}z_{2})}{[(z-z_{1})(z-z_{2})]^{n+1}}\;,$
(41)

where $\mathcal{C}^{\prime}$ is a closed anti-clockwise contour that encloses
the pole at $z=z_{1}$. but with the poles at $z=z_{0},z_{2}$ in the exterior.
For the low-priority marginal PMF, we have

$f_{\text{lo}}(n)=-\frac{(1-r)(-r_{2})^{n}}{r}\oint_{\mathcal{C}^{\prime}}\frac{dz}{2\pi
i}\,\frac{z^{n}}{z-z_{0}}{\cdot}\frac{(z^{2}-z_{1}z_{2})}{[(z-z_{1})(z-z_{2})]^{n+1}}\;.$
(42)

##  5 R-Integrals

In order to evaluate the integral representations for the joint and marginal
PMFs, derived the foregoing section, we introduce a collection of complex
contour integrals, to which we shall refer as the R-integrals, according to
the definition

$R^{m}_{n}\equiv\oint_{\mathcal{C}^{\prime}}\frac{dz}{2\pi
i}\,\frac{1}{z-z_{0}}{\cdot}\frac{z^{m}}{[(z-z_{1})(z-z_{2})]^{n}}\;,$ (43)

for $m,n=0,1,2,\ldots$, where $\mathcal{C}^{\prime}$ is a closed anti-
clockwise contour that encloses the pole at $z=z_{1}$. but with the poles at
$z=z_{0},z_{2}$ in the exterior. An immediate consequence of this definition
is the (backwards) recurrence relation

$R^{m}_{n-1}=R^{m+2}_{n}-(z_{1}+z_{2})R^{m+1}_{n}+z_{1}z_{2}R^{m}_{n}\;.$ (44)

One may also note the scaling behaviour

$R^{m}_{n}(z_{0},z_{1},z_{2})=z_{0}^{m-2n}R^{m}_{n}(1,z_{1}/z_{0},z_{2}/z_{0})\;,$
(45)

or, more generally,

$R^{m}_{n}(z_{0},z_{1},z_{2})=\zeta^{m-2n}R^{m}_{n}(z_{0}/\zeta,z_{1}/\zeta,z_{2}/\zeta)\;,$
(46)

for any $\zeta>0$.

In the present application to the priority queue, the parameters
$z_{0},z_{1},z_{2}$ are given by (39). In terms of the R-integrals, the joint
PMF is given by

$f(n,m)=\frac{(1-r)(-r_{2})^{n}}{r}\left(R^{m+n+3}_{n+1}-R^{m+n+2}_{n+1}-z_{1}z_{2}R^{m+n+1}_{n+1}+z_{1}z_{2}R^{m+n}_{n+1}\right)\;.$
(47)

If we introduce the difference functions $\Delta R^{m}_{n}\equiv
R^{m+1}_{n}-R^{m}_{n}$, then we can write

$f(n,m)=\frac{(1-r)(-r_{2})^{n}}{r}\left(\Delta
R^{m+n+2}_{n+1}-z_{1}z_{2}\Delta R^{m+n}_{n+1}\right)\;.$ (48)

Likewise, in terms of the R-integrals, we have for the low-priority marginal
PMF,

$f_{\text{lo}}(n)=-\frac{(1-r)(-r_{2})^{n}}{r}\left(R^{n+2}_{n+1}-z_{1}z_{2}R^{n}_{n+1}\right)\;,$
(49)

for $n=0,1,2,\ldots$. For the exclusively-low distribution, defined by
$f_{\text{xlo}}(n)\equiv f(n,0)$, we can write

$f_{\text{xlo}}(n)=\frac{(1-r)(-r_{2})^{n}}{r}\left(\Delta
R^{n+2}_{n+1}-z_{1}z_{2}\Delta R^{n}_{n+1}\right)\;.$ (50)

It gives the probability of finding $n$ low-priority clients in the queue and
no high-priority clients. It has a form that is similar to the low-priority
marginal $f_{\text{lo}}(n)$, and we will show later that the two are, in fact,
closely related. This relationship will provide a useful diagnostic test of
the numerical performance of the R-integral computation.

We have succeeded in recasting the problem into one that involves complex
contour integration over a collection of totally meromorphic functions. In
Figure 4, we plot the $z$-contour $\mathcal{C}^{\prime}$ that results from
taking the $p$-contour $\mathcal{C}$ to be the unit circle centred on the
origin, plotted for the case of total traffic intensity $r=0.95$ and fraction
of high-priority arrivals $\nu=0.75$. Also displayed are the locations of
R-integral poles $z_{0},z_{1},z_{2}$.

###  5.1 Recurrence

If we cast the recurrence relation (44) as

$R^{m}_{n}=R^{m-2}_{n-1}-z_{1}z_{2}R^{m-2}_{n}+(z_{1}+z_{2})R^{m-1}_{n}\;,$
(51)

for $m=2,3,\ldots$, $n=1,2,\ldots$, then it may, in principle, be solved
recursively for the $R^{m}_{n}$ starting from the seed values

$\displaystyle\begin{aligned} R^{m}_{0\hphantom{+1}}&=0\;,\\\
R^{0}_{n+1}&=\frac{(-1)^{n}}{\left[(z_{1}-z_{0})(z_{1}-z_{2})\right]^{n+1}}{\cdot}p_{n}\left(\frac{z_{1}-z_{0}}{z_{1}-z_{2}}\right)\;,\\\
R^{1}_{n+1}&=\frac{(-1)^{n}}{(z_{1}-z_{2})^{2n+1}}\binom{2n}{n}+z_{0}R^{0}_{n+1}\;,\end{aligned}$
(52)

with the polynomials $p_{n}(x)$ defined by

$p_{n}(x)\equiv\sum_{k=0}^{n}\binom{k+n}{k}x^{k}\;.$ (53)

Unfortunately, this recursion scheme is numerically unstable, especially for
small $\nu$.

###  5.2 Series Representation

Applying Cauchy’s theorem to (43), followed by an invocation of Leibniz’s
formula, we obtain

$\displaystyle\begin{aligned}
R^{m}_{n+1}&=\frac{1}{n!}{\cdot}\frac{d^{n}}{dz_{1}^{n}}\left[\frac{1}{(z_{1}-z_{2})^{n+1}}{\cdot}\frac{z_{1}^{m}}{z_{1}-z_{0}}\right]\\\
&=\sum_{k=0}^{n}\frac{1}{(n-k)!}\frac{d^{n-k}}{dz_{1}^{n-k}}\left[\frac{1}{(z_{1}-z_{2})^{n+1}}\right]{\cdot}\frac{1}{k!}\frac{d^{k}}{dz_{1}^{k}}\left[\frac{z_{1}^{m}}{z_{1}-z_{0}}\right]\;.\end{aligned}$
(54)

The first differentiation is trivial to perform, yielding

$R^{m}_{n+1}=(-1)^{n}\sum_{k=0}^{n}\binom{2n-k}{n}\frac{z_{0}^{m-k-1}}{(z_{2}-z_{1})^{2n+1-k}}S^{m}_{k}(z_{1}/z_{0})\;,$
(55)

where

$S_{k}^{m}(x)\equiv\frac{1}{k!}\frac{d^{k}}{dx^{k}}\left(\frac{x^{m}}{1-x}\right)\;.$
(56)

The functions $S_{k}^{m}(x)$ satisfy the relationship

$S_{k}^{m+1}(x)-S_{k}^{m}(x)=-\binom{m}{k}x^{m-k}\;.$ (57)

It is convenient to introduce polynomials

$P_{k}^{m}(x)\equiv(1-x)^{k+1}S_{k}^{m}(x)\;,$ (58)

so that $P_{k}^{0}(x)=1$, $P_{0}^{m}(x)=x^{m}$. Then we can write

$R^{m}_{n+1}=\frac{(-1)^{n}}{(z_{2}-z_{1})^{2n+2}}\sum_{k=0}^{n}\binom{2n-k}{n}\left(\frac{z_{2}-z_{1}}{1-z_{1}/z_{0}}\right)^{k+1}z_{0}^{m-k-1}P^{m}_{k}(z_{1}/z_{0})\;.$
(59)

Combining (56) and (58), we can establish that, for $m>k$,

$\displaystyle\begin{aligned}
P_{k}^{m}(x)&=(1-x)^{k+1}\sum_{\ell=0}^{k}\frac{1}{\ell!}\frac{d^{\ell}x^{m}}{dx^{\ell}}{\cdot}\frac{1}{(k-\ell)!}\frac{d^{k-\ell}}{dx^{k-\ell}}\left(\frac{1}{1-x}\right)\\\
&=\sum_{\ell=0}^{k}D_{\ell}^{m}(x)\;,\end{aligned}$ (60)

where

$D_{\ell}^{m}(x)\equiv\binom{m}{\ell}x^{m-\ell}(1-x)^{\ell}\;.$ (61)

Equation (61) represents a cumulative sum, each term of which can be computed
recursively. For example, when $x$ is bounded away for zero,

$D^{m}_{\ell}(x)=\left(\frac{m+1}{\ell}-1\right){\cdot}\left(\frac{1}{x}-1\right)D^{m}_{\ell-1}(x)\;,$
(62)

for $\ell=1,2,\ldots$, with $D^{m}_{0}(x)=x^{m}$. A similar recursion holds
for small $x$, computed backwards from $D^{m}_{m}(x)=(1-x)^{m}$.

An explicit representation of the polynomials $P_{k}^{m}(x)$ is given by

$P_{k}^{m}(x)=1-(1-x)^{k+1}\sum_{\ell=0}^{m-k-1}\binom{k+\ell}{\ell}x^{\ell}\;.$
(63)

It may be observed that $P_{k}^{m}(x)=1$ whenever $m\leq k$, and that
$P_{k}^{m}(x)\geq 0$ for all $0\leq x\leq 1$. These polynomials also satisfy
the recurrence relation

$P_{k}^{m+1}(x)=P_{k}^{m}(x)+\frac{m}{k}(1-x)\left[P_{k-1}^{m}(x)-P_{k-1}^{m-1}(x)\right]\,$
(64)

for $k,m=1,2,\ldots$, subject to

$P_{0}^{m}(x)=x^{m}\;,\quad P_{k}^{0}(x)=1\;,\quad
P_{k}^{1}(x)=1-(1-x)\delta_{k0}\;.$ (65)

###  5.3 Evaluation

In order to achieve good numerical behaviour as $\nu\to 1$, it is convenient
to work with the scaled integrals
$\hat{R}^{m}_{n+1}\equiv(-r_{\text{lo}})^{n}R^{m}_{n+1}$, for which we have
the well-behaved series representation

$\begin{split}\hat{R}^{m}_{n+1}=&\frac{1}{(z_{2}-z_{1})(1-z_{1}/z_{0})}{\cdot}\left(\frac{r_{\text{lo}}}{(z_{2}-z_{1})^{2}}\right)^{n}\\\
&{}\times\sum_{k=0}^{n}\binom{2n-k}{n}\left(\frac{z_{2}-z_{1}}{1-z_{1}/z_{0}}\right)^{k}z_{0}^{m-k-1}P_{k}^{m}(z_{1}/z_{0})\;.\end{split}$
(66)

Thus, we consider the computation of the vectors

$\hat{\mathbf{R}}^{(m)}\equiv[\hat{R}^{m}_{1},\hat{R}^{m}_{2},\ldots,\hat{R}^{m}_{N+1}]^{\sf
T}\;.$ (67)

To assist with this, we define the constant

$\kappa\equiv\frac{1}{(z_{2}-z_{1})(1-z_{1}/z_{0})}\;,$ (68)

the diagonal matrices

$\displaystyle\begin{aligned} A&\equiv\mathop{\rm
diag}\nolimits[a^{0},a^{1},\ldots,a^{N}]\;,&a&\equiv
r_{\text{lo}}/(z_{2}-z_{1})^{2}\;,\\\ B&\equiv\mathop{\rm
diag}\nolimits[b^{0},b^{1},\ldots,b^{N}]\;,&b&\equiv(z_{2}-z_{1})/(1-z_{1}/z_{0})\;,\end{aligned}$
(69)

and the combinatorial matrix

$C_{nk}\equiv\binom{2n-k}{n}$ (70)

provided $k\leq n$ and is zero otherwise. We also introduce the polynomial
vectors

$\mathbf{P}^{(m)}\equiv[P^{(m)}_{0},P^{(m)}_{1},\ldots,P^{(m)}_{N}]^{\sf
T}\;,\quad P^{(m)}_{k}\equiv z_{0}^{m-k-1}P_{k}^{m}(z_{1}/z_{0})\;.$ (71)

Then, we can write (66) as

$\displaystyle\begin{aligned}
\hat{\mathbf{R}}^{(m)}&=\kappa{\cdot}ACB\mathbf{P}^{(m)}\\\
&=\kappa{\cdot}(ACA^{-1}){\cdot}AB{\cdot}\mathbf{P}^{(m)}\;.\end{aligned}$
(72)

At this point, we note that the product $AB$ is the diagonal matrix of
increasing powers

$AB=\mathop{\rm
diag}\nolimits[\gamma^{0},\gamma^{1},\ldots,\gamma^{N}]\;,\quad\gamma\equiv
r_{\text{lo}}/[(z_{2}-z_{1})(1-z_{1}/z_{0})]\;,$ (73)

and that $(ACA^{-1})_{nk}=a^{n-k}C_{nk}$, which is easily computed by
observing the cumulative product form

$a^{\ell}\binom{n+\ell}{n}=\prod_{j=1}^{\ell}\left(1+n/j\right)a\;.$ (74)

If we combine the column vectors $\hat{\mathbf{R}}^{(m)}$ and
$\mathbf{P}^{(m)}$ into respective matrices, so that

$\displaystyle\begin{aligned}
\hat{\mathbf{R}}&\equiv[\hat{\mathbf{R}}^{(0)},\hat{\mathbf{R}}^{(1)},\ldots,\hat{\mathbf{R}}^{(M)}]\;,\\\
\mathbf{P}&\equiv[\mathbf{P}^{(0)},\mathbf{P}^{(1)},\ldots,\mathbf{P}^{(M)}]\;,\end{aligned}$
(75)

then we obtain the matrix equation

$\hat{\mathbf{R}}=\kappa{\cdot}(ACA^{-1}){\cdot}AB{\cdot}\mathbf{P}\;.$ (76)

In Figure 5.4, we plot the queue-length PMF for the low-priority arrivals, as
the negative base-10 logarithm, for total traffic intensity $r=0.99$ and a
range of hifrac values $\nu$. Overlaid, are the asymptotic curves in the large
queue-length limit. This is given by

$f_{\text{lo}}(n)\mathrel{\mathop{\sim}_{n\to\infty}}\sqrt{\frac{1-r}{\pi
r}}{\cdot}\frac{r^{n}}{\sqrt{n}}\;,$ (77)

when $r_{\text{hi}}=r^{2}$ (or equivalently $\nu=r$). Otherwise, the low-
priority marginal PMF can be decomposed into two components according to

$f_{\text{lo}}(n)=f_{\text{pol}}(n){\cdot}\Theta(r^{2}-r_{\text{hi}})+f_{\text{cut}}(n)\;,$
(78)

where $\Theta(x)$ denotes the Heaviside function such that $\Theta(x)=1$ for
$x\geq 0$ and vanishes otherwise. The large-$n$ behaviour of these components
is given by

$\displaystyle\begin{aligned}
f_{\text{pol}}(n)&\mathrel{\mathop{\sim}_{n\to\infty}}\left[1-\frac{r(1-r)}{r_{\text{lo}}}\right](1-r)r^{n-1}\;,\\\
f_{\text{cut}}(n)&\mathrel{\mathop{\sim}_{n\to\infty}}\frac{(\sqrt{r_{\text{hi}}}/r_{\text{lo}})^{1/2}}{2\sqrt{\pi}r}{\cdot}\frac{1-r}{(\chi-1/r)\chi^{n-1/2}n^{3/2}}\;,\end{aligned}$
(79)

where $\chi\equiv 1+(1-\sqrt{r_{\text{hi}}})^{2}/r_{\text{lo}}>1/r$. The
derivation of these results, which will be presented in a forthcoming paper,
follows directly from the pole/cut integral representation of the
distribution, mentioned in Section 4. The computed points, represented by the
coloured dots, are interpolated by black curves. The asymptotic curves are
indicated by a coloured dashed line-style. Thus, when the interpolation
between the data points becomes coloured, this indicates that the agreement
between the computation and asymptotic limit is within the linewidth of the
graph. In Figure 5.4, we plot the queue-length PMF for the low-priority
arrivals, as the negative base-10 logarithm, for the case of total traffic
intensity $r=0.99$ and fraction of high-priority arrivals $\nu=0.95$, where
asymptotic behaviour is slow to set in. We see that the computation remains
robust up to a queue length of at least $n=1000$ which lies deep in the
asymptotic region. In Figure 5.4, we plot a two-dimensional map of the joint
probability distribution $f(n,m)$ of the queue lengths, for total traffic
intensity $r=0.75$ and fraction of high-priority arrivals $\nu=0.9$. A
logarithmic scaling has been applied, such that
$f(n,m)\leftarrow\max\\{0,1+\log_{10}(f(n,m)/f_{\text{max}})/20\\}$, where
$f_{\text{max}}\equiv\max\\{f(n,m)\\}$.

###  5.4 Limiting Cases

The $\nu\to 0$ limiting behaviour of the R-integrals is given by

$R_{n}^{m}\mathrel{\mathop{\sim}_{\nu\to
0^{+}}}\left\\{\begin{array}[]{cll}0&\quad\text{for}&m>n\\\
(-1)^{n-1}\left[1-r^{n-1}/(1+r)^{2n-1}\right]&\quad\text{for}&m=n\\\
(-1)^{n-1}/\nu^{n-m}&\quad\text{for}&m<n\end{array}\right.\;.$ (80)

At the opposite extreme, for $\nu=1$, we have $z_{0}=1,\;z_{1}=r,\;z_{2}=1$,
in which case

$\frac{z_{2}-z_{1}}{1-z_{1}/z_{0}}=1\;.$ (81)

It follows that

$R_{n+1}^{m}\mathrel{\mathop{\sim}_{\nu\to
1^{-}}}\frac{(-1)^{n}}{(1-r)^{2n+2}}\sum_{k=0}^{n}\binom{2n-k}{n}P_{k}^{m}(r)\;.$
(82)

Equation (80) shows that the R-integrals become singular for small $\nu$ when
$m<n$. This is one reason for the numerical instability of the recurrence
relations (51), given that the seed values always reside in this region.

x Queue-length PMF for the low-priority arrivals, plotted as the negative
base-10 logarithm, for total traffic intensity $r=0.9$ and a range of hifrac
values ($\nu$). Asymptotic curves for the large queue-length limit are
overlaid.

x Queue-length PMF for the low-priority arrivals, plotted as the negative
base-10 logarithm, for total traffic intensity $r=0.99$ and fraction of high-
priority arrivals $\nu=0.9$, with queue lengths extending far into the
asymptotic region. It is compared with the exact asymptotic curve in the large
queue-length limit.

x 2D map of the joint probability distribution of the queue lengths, with a
logarithmic scaling applied, for total traffic intensity $r=0.75$ and fraction
of high-priority arrivals $\nu=0.9$.

##  6 Numerical Tests

x Four tests of the joint probability distribution of the queue lengths as
functions of the fraction of high-priority arrivals $\nu$, across a wide range
of values for the total traffic intensity $r$, as displayed. The values on the
vertical axes indicate the number of decimal places of agreement.

Various tests can be applied to quantify the numerical performance of the
algorithm for the computation of the joint PMF.

###  6.1 Aggregation Test

The aggregated queue-length distribution describes the total number of
entities in the queue, regardless of priority level. This is equivalent to the
queue-length distribution of the basic M/M/$c$ queueing model with traffic
intensity $r=r_{\text{lo}}+r_{\text{hi}}$, which is known to be a simple
geometric distribution. Hence, the exact aggregate PMF is given by

$f_{\text{agg}}^{(\text{ex})}(k)=(1-r)r^{k}\;,$ (83)

for $k=0,1,2,\ldots$

One diagnostic test of the R-integral computational methodology is to check
how well the aggregate PMF constructed from the computed joint PMF reproduces
the exact result. This test is more convenient than similarly testing against
the marginals as only a finite summation is required. Considering the joint
PMF as a matrix whose rows and columns are labelled by its integer arguments,
values of the aggregate PMF are given by successive finite sums along the
anti-diagonals. Specifically, in terms of the R-integrals, the aggregate PMF
is expressed as

$\displaystyle\begin{aligned} f_{\text{agg}}(k)&=\sum_{n=0}^{k}f(n,k-n)\\\
&=\frac{1-r}{r}\sum_{n=0}^{k}(-r_{\text{lo}})^{n}(\Delta
R_{n+1}^{k+2}-r_{\text{hi}}\Delta R_{n+1}^{k})\;,\end{aligned}$ (84)

for $k=0,1,2,\ldots$.

We then consider the measure of performance (MOP)

$\Xi_{\text{agg}}\equiv-\max_{k\geq
0}\left\\{\log_{10}\left(|\ln(f_{\text{agg}}(k))-\ln(f_{\text{agg}}^{(\text{ex})}(k))|\right)\right\\}\;,$
(85)

where the maximum is taken over all values $0\leq k\leq n_{\text{lim}}$ such
that $f_{\text{agg}}^{(\text{ex})}(k)>p_{\text{lim}}>0$. Since we are working
in double-precision arithmetic444All computation is performed in Matlab
R2020a, which implements IEEE Standard 754 for double precision., all MOPs of
this kind are capped at a maximum allowed value of $16$. The interpretation of
$\Xi_{\text{agg}}$ (and similarly for all of the subsequent MOPs) is that it
indicates the number of decimal places of numerical agreement in the worst
case.

###  6.2 Xhi-Test

The exclusively-high distribution, defined by $f_{\text{xhi}}(m)\equiv
f(0,m)$, gives the probability of finding $m$ high-priority clients in the
queue and no low-priority clients. An exact expression for the exclusively-
high probability is given by

$f_{\text{xhi}}^{(\text{ex})}(m)=(1-r)(r_{\text{hi}}/z_{2})^{m}\;,$ (86)

for $m=0,1,2,\ldots$ It is simple to calculate directly as the R-integral has
only a simple pole when $n=0$. One should note that $f_{\text{xhi}}(m)$ is not
a proper PMF since $\sum_{m=0}^{\infty}f_{\text{xhi}}(m)<1$, unless $\nu=1$,
but can be turned into a conditional PMF by means of an overall scale factor.

In terms of the R-integrals, the exclusively-high PMF is expressed as

$f_{\text{xhi}}(m)=\frac{1-r}{r}\left(\Delta R^{m+2}_{1}-z_{1}z_{2}\Delta
R^{m}_{1}\right)\;,$ (87)

and we consider the MOP

$\Xi_{\text{xhi}}\equiv-\max_{m\geq
0}\left\\{\log_{10}\left(|\ln(f_{\text{xhi}}(m))-\ln{(}f_{\text{xhi}}^{(\text{ex})}(m))|\right)\right\\}\;,$
(88)

where the maximum is taken over all values $0\leq m\leq n_{\text{lim}}$ such
that $f_{\text{xhi}}^{(\text{ex})}(m)>p_{\text{lim}}>0$.

###  6.3 Xlo-Test

Checking whether the computed joint PMF gives rise to the correct marginal
distribution, numerically, is not a convenient enterprise as it necessitates
an infinite summation. However, it is possible to devise an alternative test
that checks the consistency of the numerical low-priority marginal with the
numerically computed joint PMF. In the xlo-test, we relate the exclusively-low
distribution $f_{\text{xlo}}(n)$ with the low priority marginal
$f_{\text{lo}}(n)$. To achieve this, we consider the PGF (8) recast into the
form

$g_{m}(p)=(1-r)\frac{1-\lambda_{1}(p)}{\lambda_{2}(p)-r}{\cdot}\lambda_{1}^{m}(p)\;.$
(89)

Specialized to the case $m=0$, this may be expressed as

$g_{0}(p)=(1-r)\left[1+r_{\text{lo}}p{\cdot}\frac{1}{\lambda_{2}(p)-r}\right]\;.$
(90)

Since the PGF of the low-priority marginal is given by

$g_{\text{lo}}(p)=\sum_{m=0}^{\infty}g_{m}(p)=\frac{1-r}{\lambda_{2}(p)-r}\;,$
(91)

we arrive at the result

$g_{0}(p)=1-r+r_{\text{lo}}p{\cdot}g_{\text{lo}}(p)\;.$ (92)

There is a generalization of this result to non-zero values of $m$ that
relates $g_{m}(p)$ to $g_{\text{lo}}(p)$. Its derivation is presented in the
Appendix. From the relationships

$g_{0}(p)=\sum_{n=0}^{\infty}p^{n}f_{\text{xlo}}(n)\;,\quad
g_{\text{lo}}(p)=\sum_{n=0}^{\infty}p^{n}f_{\text{lo}}(n)\;,$ (93)

we can equate powers to read off that

$f_{\text{xlo}}(0)=1-r\;,\quad
f_{\text{xlo}}(n)=r_{\text{lo}}{\cdot}f_{\text{lo}}(n-1)\;,$ (94)

for $n=1,2,\ldots$, or, equivalently,

$f_{\text{xlo}}(n)=(1-r)\delta_{n0}+(1-\delta_{n0})r_{\text{lo}}{\cdot}f_{\text{lo}}(n-1)\;,$
(95)

for $n=0,1,2,\ldots$, where we can formally set $f_{\text{lo}}(-1)\equiv 0$.
One should note that $f_{\text{xlo}}(n)$ is not a proper PMF since

$\sum_{n=0}^{\infty}f_{\text{xlo}}(n)=1-r_{\text{hi}}\;,$ (96)

but can be turned into a conditional PMF by means of an overall scale factor.

For the xlo-test, the LHS of (95) is taken to be given by (50) and is compared
with the RHS of (95) where the marginal PMF $f_{\text{lo}}(n)$ is expressed in
terms of the R-integrals via (49). The relevant MOP is taken to be

$\Xi_{\text{xlo}}\equiv-\max_{n>0}\left\\{\log_{10}\left(|\ln(f_{\text{xlo}}(n))-\ln{(}r_{\text{lo}}f_{\text{lo}}(n-1))|\right)\right\\}\;,$
(97)

where the maximum is taken over all values $0<n\leq n_{\text{lim}}$ such that
$f_{\text{xlo}}(m)>p_{\text{lim}}>0$.

###  6.4 Nearest-Neighbour Test

A direct consequence of the recurrence relations for the R-integrals is that
the joint PMF at any given interior point $(n,m)$ is a positively weighted sum
of the joint PMF values at three of its four nearest neighbours:

$f(n,m)=\frac{1}{1+r}\left[f(n,m+1)+r_{\text{lo}}f(n-1,m)+r_{\text{hi}}f(n,m-1)\right]\;,$
(98)

for all $m,n>0$.

In order to apply the neighbour test, we first compute the joint PMF $f(n,m)$
on a 2D grid of points $(n,m)$ from (48). Next, we use these values to compute
the RHS of (98), which we shall denote $f_{\text{nn}}(n,m)$. Then, we consider
the MOP

$\Xi_{\text{nn}}\equiv-\max_{m,n>0}\left\\{\log_{10}\left(|\ln(f(n,m))-\ln(f_{\text{nn}}(n,m))|\right)\right\\}\;,$
(99)

where the maximum is taken over all values $0<m,n\leq n_{\text{lim}}$ such
that $f(n,m)>p_{\text{lim}}>0$.

###  6.5 Quadratic Test

In this test, we compare the results for the joint queue-length PMF computed
from the R-integral (denoted $f_{\text{ri}}(n,m)$) with that computed by the
quadratic recurrence (denoted $f_{\text{qr}}(n,m)$). The MOP is taken to the
be number of decimal places of agreement, as given by

$\Xi_{\text{qr}}\equiv-\max_{m,n>0}\left\\{\log_{10}\left(|\ln(f_{\text{ri}}(n,m))-\ln(f_{\text{qr}}(n,m))|\right)\right\\}\;,$
(100)

where the maximum is taken over all values $0<m,n\leq n_{\text{lim}}$ such
that $f(n,m)>p_{\text{lim}}>0$.

###  6.6 Results

Figure 6 presents the results of the numerical tests. The MOP values relevant
to the R-integral computations are displayed on the vertical axis against the
full range of high-priority arrival fraction (hifrac) $\nu$ on the horizontal
axis. Individual curves are plotted for a discrete collection of traffic
intensities, spanning a wide range. Agreement always exceeds eight decimal
places, and is generally much higher. The nearest-neighbour and xlo-tests
check the internal consistency of the computations, while the aggregation and
xhi-tests check against exact analytical results.

The maximum queue occupancy to be examined was taken to be
$n_{\text{lim}}=1000$. PMF intervals examined included everything down to a
tail value of $p_{\text{lim}}=10^{-20}$ except in the xhi-test where
$p_{\text{lim}}=10^{-30}$ was used.

Figure 6.6 presents the results of comparing the joint queue-length
distribution computed from the R-integral with that computed by the quadratic
recurrence. The close agreement observed implies a high level of accuracy for
each method across the complete range of parameters. Worst case accuracy
occurs when both the traffic intensity $r$ and hifrac $\nu$ approach unity. In
Table 6.6, we present results that investigate this region in more detail.
Values of $r$ close to unity have been reported to be problematic for the
matrix-geometric approach [17, 18]. The table shows that both of the present
methods behave well in this region. The fourth and fifth columns indicate the
smallest rectangular subset $[0,n_{\text{lo}}]\times[0,n_{\text{hi}}]$ of
$[0,n_{\text{lim}}]\times[0,n_{\text{lim}}]$ that contains all grid points
$(m,n)$ with probability greater than $p_{\text{lim}}=10^{-20}$. A value of
$n_{\text{lim}}=1000$ in one or both columns indicates that $p_{\text{lim}}$
was not attained in some direction. The last column is the minimum probability
that was achieved over all considered grid points whose probability values
exceed $p_{\text{lim}}$. Computation time for the quadratic recurrence method
is two orders of magnitude faster than for the R-integral method. Finally,
Figure 6.6 repeats the quadratic test as described above, but for the low-
priority marginal PMFs, with the distribution arising from the quadratic
recurrence computed by the algorithm of (23). The legend indicates the of
range maximum queue lengths $n$ that had to be considered across the full
range of hifrac values $\nu$ in order to the attain the limiting probability
level $p_{\text{lim}}=10^{-20}$ for the given traffic intensity $r$. Agreement
between the R-integral and quadratic recurrence approaches is observed to
exceed ten decimal places in the worst case..

Joint PMF Comparison  $r$ $\nu$ $\Xi_{\text{qr}}$ $n_{\text{hi}}$
$n_{\text{lo}}$ $p_{\text{min}}$ $0.99$ $0.95$ $9.3279$ $609$ $1000$
$1.0000\times 10^{-20}$ $0.99$ $8.1611$ $1000$ $1000$ $1.0000\times 10^{-20}$
$0.999$ $6.6633$ $1000$ $1000$ $1.0000\times 10^{-20}$ $1.00$ $11.7428$ $1000$
$0$ $4.3171\times 10^{-7}$ $0.999$ $0.95$ $9.4247$ $685$ $1000$ $1.0000\times
10^{-20}$ $0.99$ $8.4169$ $1000$ $1000$ $1.0017\times 10^{-20}$ $0.999$
$7.2251$ $1000$ $1000$ $6.6926\times 10^{-18}$ $1.00$ $9.6972$ $1000$ $0$
$3.6770\times 10^{-4}$ $0.9999$ $0.95$ $9.4344$ $657$ $1000$ $1.0000\times
10^{-20}$ $0.99$ $8.4361$ $1000$ $1000$ $1.0000\times 10^{-20}$ $0.999$
$7.2455$ $1000$ $1000$ $1.0540\times 10^{-18}$ $1.00$ $7.8504$ $1000$ $0$
$9.0483\times 10^{-5}$

x Comparison of the joint probability distribution of the queue lengths as
computed from the R-integral and from the quadratic recurrence. The MOP is
plotted as a function of the fraction of high-priority arrivals $\nu$, across
a wide range of values for the total traffic intensity $r$, as displayed. The
values on the vertical axes indicate the number of decimal places of
agreement.

x Comparison of the low-priority marginal distribution of the queue lengths as
computed from the R-integral and from the quadratic recurrence. The MOP is
plotted as a function of the fraction of high-priority arrivals $\nu$, across
a wide range of values for the total traffic intensity $r$, as displayed. The
values on the vertical axes indicate the number of decimal places of
agreement.

The exact results for the queue-lengths distributions derived here were also
tested against Monte-Carlo simulation. Excellent agreement was found across
the entire parametric domain. Details will be presented elsewhere. We have
also checked against the results in Table 3 of [8] (where the service times
are equal) to find complete agreement. There, in the case of the present
problem, the quantity $P_{\text{Q}}$ is related to the no-wait probability
$P_{\text{NW}}$ given in (3) by $P_{\text{Q}}=1-P_{\text{NW}}$, and $p(0,0)$
is the probability that the system is empty, given by

$\frac{1}{p(0,0)}=\frac{(rN)^{N}}{N!}\left[\frac{1}{1-r}+\Gamma_{\text{scl}}(rN,N)\right]\;,$
(101)

where

$\Gamma_{\text{scl}}(x,\nu)\equiv\frac{\nu e^{x}}{x^{\nu}}\int_{x}^{\infty}dt\
t^{\nu-1}e^{-t}$ (102)

is the scaled upper incomplete gamma function as implemented in Matlab. In the
present problem, neither of these quantities depend on the priority structure.
We relate the mean waiting times given in the table to the mean queue lengths
via Little’s law.

##  7 Conclusions

Simple methods for highly accurate computation of the joint and marginal
distributions for a non-preemptive two-level priority queue have been
developed. Explicit closed-form representations for the joint and marginal
PMFs have also been derived, something that has not been achieved previously.
Future work could entail extension of the present methods to unequal services
rates among the priority levels.

Convolutional Form In this appendix, we derive a relationship between joint
PMF $f(n,m)$ and the low-priority marginal $f_{\text{lo}}(n)$. One may observe
the general structure

$g_{m}(p)=A^{(m)}(p)+B^{(m)}(p){\cdot}g_{\text{lo}}(p)\;,$ (103)

for polynomials

$A^{(m)}(p)=\sum_{n=0}^{m}A^{(m)}_{n}p^{n}\;,\quad
B^{(m)}(p)=\sum_{n=0}^{m}B^{(m)}_{n}p^{n}\;.$ (104)

Thus, we obtain the convolutional form for the joint PMF:

$f(n,m)=A^{(m)}_{n}+\sum_{k=0}^{n}B^{(m)}_{k}{\cdot}f_{\text{lo}}(n-k)\;,$
(105)

which generalizes the relationship between $f_{\text{xlo}}$ and
$f_{\text{lo}}$ given in (95) to non-zero values of $m$. One may also note the
special case

$f_{\text{xhi}}(m)\equiv f(0,m)=A^{(m)}_{0}+B^{(m)}_{0}f_{\text{lo}}(0)\;.$
(106)

In what follows, we derive explicit expressions for $A^{(m)}(p)$ and
$B^{(m)}(p)$.

Since $\lambda_{1}(p)$ satisfies a quadratic equation, we have that
$\lambda_{1}^{m}(p)=\alpha_{m}(p)\lambda_{1}(p)+\beta_{m}(p)$ for some
polynomials $\alpha_{m}(p)$, $\beta_{m}(p)$. The fact that
$\alpha_{m},\beta_{m}$ are polynomials follows from examining
$\lambda_{1}^{2}$. Let us now recall that the Chebyshev polynomials of the
first and second kind, $T_{n}(x)$ and $U_{n}(x)$ respectively, may be
expressed as

$\displaystyle\begin{aligned}
T_{n}(x)&=\tfrac{1}{2}\left[(x+\sqrt{x^{2}-1})^{n}+(x-\sqrt{x^{2}-1})^{n}\right]\;,\\\
\sqrt{x^{2}-1}U_{n-1}(x)&=\tfrac{1}{2}\left[(x+\sqrt{x^{2}-1})^{n}-(x-\sqrt{x^{2}-1})^{n}\right]\;,\end{aligned}$
(107)

which implies the identity

$(x-\sqrt{x^{2}-1})^{n}=T_{n}(x)-xU_{n-1}(x)+(x-\sqrt{x^{2}-1}){\cdot}U_{n-1}(x)\;.$
(108)

Since we have

$\lambda_{1}(p)=b(p)-\sqrt{b^{2}(p)-r_{1}}\;,\quad
b(p)\equiv(1+r-r_{2}p)/2\;,$ (109)

it follows that

$\frac{\lambda_{1}^{m}(p)}{r_{1}^{m/2}}=T_{m}(x(p))-x(p){\cdot}U_{m-1}(x(p))+\frac{\lambda_{1}(p)}{\sqrt{r_{1}}}{\cdot}U_{m-1}(x(p))\;,$
(110)

with $x(p)\equiv b(p)/\sqrt{r_{1}}$. Consequently,

$\displaystyle\begin{aligned}
\alpha_{m}(p)&=r_{1}^{(m-1)/2}U_{m-1}(x(p))\;,\\\
\beta_{m}(p)&=r_{1}^{m/2}\left[T_{m}(x(p))-x(p){\cdot}U_{m-1}(x(p))\right]\;,\end{aligned}$
(111)

for $x=0,1,\ldots$, with $U_{-1}(x)\equiv 0$. The first few
$\alpha$-coefficients are given by

$\alpha_{0}(p)=0\;,\quad\alpha_{1}(p)=1\;,\quad\alpha_{2}(p)=1+r-r_{2}p\;.$
(112)

The first few $\beta$-coefficients are given by

$\beta_{0}(p)=1\;,\quad\beta_{1}(p)=0\;,\quad\beta_{2}(p)=-r_{1}\;.$ (113)

Using identities satisfied by the Chebyshev polynomials, (111) can be
simpified as

$\alpha_{m}(p)=r_{1}^{(m-1)/2}U_{m-1}(x(p))\;,\quad\beta_{m}(p)=-r_{1}^{m/2}U_{m-2}(x(p))\;,$
(114)

which implies the relationship $\beta_{m}(p)=-r_{1}\alpha_{m-1}(p)$, where we
formally set $\alpha_{-1}(p)\equiv 1$.

Noting that

$\displaystyle\begin{aligned}
(1-\lambda_{1})\lambda_{1}^{m}&=(\alpha_{m}-\alpha_{m+1})\lambda_{1}+(\beta_{m}-\beta_{m+1})\\\
&=-(\alpha_{m}-\alpha_{m+1})(\lambda_{2}-r)+(1-r_{2}p)(\alpha_{m}-\alpha_{m+1})+(\beta_{m}-\beta_{m+1})\;,\end{aligned}$
(115)

followed by substitution into the representation

$g_{m}(p)=\frac{1-r}{\lambda_{2}-r}{\cdot}(1-\lambda_{1})\lambda_{1}^{m}\;,$
(116)

yields

$g_{m}(p)=-(1-r)(\alpha_{m}-\alpha_{m+1})+g_{\text{lo}}(p)\left[(1-r_{2}p)(\alpha_{m}-\alpha_{m+1})+(\beta_{m}-\beta_{m+1})\right]\;,$
(117)

from which we can read off

$\displaystyle\begin{aligned}
A^{(m)}(p)&=(1-r)\left[\alpha_{m+1}(p)-\alpha_{m}(p)\right]\;,\\\
B^{(m)}(p)&=(1-r_{2}p)\left[\alpha_{m}(p)-\alpha_{m+1}(p)\right]+\left[\beta_{m}(p)-\beta_{m+1}(p)\right]\;.\end{aligned}$
(118)

In terms of Chebyshev polynomials, this becomes

$\displaystyle\begin{aligned}
A^{(m)}(p)&=-(1-r)r_{1}^{(m-1)/2}\left[U_{m-1}(x(p))-\sqrt{r_{1}}U_{m}(x(p))\right]\;,\\\
B^{(m)}(p)&=-\frac{1-r_{2}p}{1-r}{\cdot}A^{(m)}(p)+\frac{r_{1}}{1-r}{\cdot}A^{(m-1)}(p)\;,\end{aligned}$
(119)

for $m=0,1,\ldots$, where we set $U_{-1}(x)\equiv 0$, $U_{-2}(x)\equiv-1$. As
a sanity check, it is straightforward to confirm that
$A^{(m)}(1)+B^{(m)}(1)=(1-r_{1})r_{1}^{m}$. Another check is given by

$\displaystyle\begin{aligned}
\sum_{m=0}^{\infty}A^{(m)}(p)&=\left\\{\begin{array}[]{cll}(1-r)/(1-r_{1})&\quad\text{for}&p=1\\\
0&\quad\text{for}&p\neq 1\end{array}\right.\;,\\\
\sum_{m=0}^{\infty}B^{(m)}(p)&=\left\\{\begin{array}[]{cll}(r-r_{1})/(1-r_{1})&\quad\text{for}&p=1\\\
1&\quad\text{for}&p\neq 1\end{array}\right.\;.\\\ \end{aligned}$ (120)

## Acknowledgments.

The authors gratefully acknowledge useful discussions with Dr. Stephen
Bocquet.

## References

  * Bertsimas and Nakazato [1995] Bertsimas D, Nakazato D (1995) The distributional Little’s law and its applications. _Operations Research_ 43(2):298–310.
  * Cidon and Sidi [1990] Cidon I, Sidi M (1990) Recursive computation of steady-state probabilities in priority queues. _Operations Research Letters_ 9(4):249–256.
  * Cobham [1954] Cobham A (1954) Priority assignment in waiting-line problems. _Operations Research_ 2(1):70–76.
  * Cobham [1955] Cobham A (1955) Priority assignment – A correction. _Operations Research_ 3(4):547–547.
  * Cohen [1956] Cohen J (1956) Certain delay problems for a full availability trunk group loaded by two traffic sources. _Philips Telecommunications Review_ 16(3):105–113.
  * Davis [1966] Davis R (1966) Waiting time distribution of a multi-server, priority queuing sytem. _Operations Research_ 14(1):133–136.
  * Dressin and Reich [1957] Dressin S, Reich E (1957) Priority assignment in a waiting line. _Quarterly of Applied Mathematics_ 15(2):208–211.
  * Gail et al. [1988] Gail H, Hantler S, Taylor B (1988) Analysis of non-preemptive priority multi-server queue. _Advances in Applied Probability_ 20(4):852–879.
  * Holley [1954] Holley J (1954) Waiting line subject to priorities. _Operations Research_ 2(3):341–343.
  * Hou and Zhao [2020] Hou J, Zhao X (2020) Using a priority queueing approach to improve emergency department performance. _Journal of Management Analytics_ 7(1):28–43.
  * Kao and Narayanan [1990] Kao E, Narayanan K (1990) Computing steady state probabilities of a non-preemptive priority queue. _INFORMS Journal on Computing_ 2(3):211–218.
  * Kao and Wilson [1999] Kao E, Wilson S (1999) Analysis of non-preemptive priority queues with multpile servers and two priority classes. _European Journal of Operational Research_ 118(1):181–193.
  * Keilson and Servi [1988] Keilson J, Servi L (1988) A distributional form of Little’s law. _Operations Research Letters_ 7(5):223–227.
  * Kella and Yechiali [1985] Kella O, Yechiali U (1985) Waiting times in the non-preemptive priority M/M/$c$ queue. _Communications in Statistics. Stochastic Models_ 1(2):257–262.
  * Kesten and Runnenberg [1957] Kesten H, Runnenberg T (1957) Priority in waiting line problems. _Koninkl. Nederl. Akademie van Wetenschappen – Amsterdam_ A60(3):312–336.
  * Marks [1973] Marks B (1973) State probabilities of M/M/$1$ priority queues. _Operations Research_ 21(4):974–987.
  * Miller [1981] Miller D (1981) Computation of steady-state probabilities for M/M/$1$ priority queues. _Operations Research_ 29(5):945–958.
  * Miller [1982] Miller D (1982) Steady-state algorithmic analysis of M/M/$c$ two-priority queues with heterogeneous rates. Disney R, Ott T, eds., _Applied Probability – Computer Science: The Interface_ , volume 3 of _Progress in Computer Science_ , 207–222 (New York, NY, USA: Springer Science + Business Media, LLC).
  * Miller [1960] Miller R (1960) Priority queues. _Annals of Mathematical Statistics_ 31(1):86–103.
  * Neuts [1981] Neuts M (1981) _Matrix-Geometric Solutions in Stochastic Models_ (Baltimore, MD, USA: Johns Hopkins University Press).
  * Orman [1995] Orman A (1995) _Models for Scheduling a Multifunction Phased Array Radar System_. Ph.D. thesis, University of Southampton, Southampton, UK.
  * Orman et al. [1996] Orman A, Potts C, Shahani A, Moore A (1996) Scheduling for a multifunction phased array radar system. _European Journal of Operational Research_ 90(1):13–25.
  * Pestalozzi [1964] Pestalozzi G (1964) Priority rules for runway use. _Operations Research_ 12(6):941–950.
  * Shortle et al. [2018] Shortle J, Thompson J, Gross D, Harris C (2018) _Fundamentals of Queueing Theory_ (Hoboken, NJ, USA: John Wiley and Sons, Inc.), fifth edition.
  * Taylor and Templeton [1980] Taylor I, Templeton J (1980) Waiting time in a multi-server cutoff-priority queue, and its application to an urban ambulance service. _Operations Research_ 28(5):1168–1188.
  * Wagner [1997] Wagner D (1997) Waiting times for a finite-capacity multi-server model with non-preemptive priorities. _European Journal of Operational Research_ 102(1):227–241.