0

Zuk and Kirszenblat

Non-Preemptive Multi-Level Priority Queue

Joint Queue-Length Distribution for the Non-Preemptive Multi-Server Multi-
Level Markovian 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 obtained using simple and exact methods for the joint
queue-length distribution of the M/M/$c$ queue with an arbitrary number of
non-preemptive priority levels. This work is the first to provide explicit
results for the joint probability generating function and joint probability
mass function for a general number of priority levels. A fixed-point iteration
is developed for the stationary balance equations, which enables direct
computation of the joint queue-length distribution. A multi-variate
probability generating function is also derived, from which the joint
probability mass function can be computed by means of a multi-dimensional fast
Fourier transform method.

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

Date created: July 27, 2023. Last update: October 24, 2023.

##  1 Introduction

This work is concerned with the development of practical algorithms for
computing the joint queue-length distribution for the non-preemptive Markovian
priority queue with a general number of priority levels. In the most recent
edition of their textbook, Shortle et al. [11] remark 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’. Elmelegy [6] has also commented that ‘given the
immense literature studying non-preemptive priority queuing systems, it is
hard to find a simple and exact method that calculates the performance
measures of non-preemptive priority systems with more than two priority
levels’. The present discussion serves to fill this knowledge gap.

Previous work on the non-preemptive priority M/M/$c$ queue has been reviewed
recently in [14]. The vast majority of effort concerning the joint
distribution has focused on the two-level problem. The single source of
previous work addressing the joint queue-length distribution for more than two
priority levels comprises the papers of Wignall [12] and Wignall and Enns
[13]. These consider single-server systems involving multiple queues with
distinct arrival rates ranked by priority level, and with feedback – requiring
deterministic or probabilistic transitions between queues before system exit
is achieved. The present problem can, in principle, be constructed as a
special case of this scheme. In [12, 13], equations that must be solved
recursively are presented for the multivariate probability generating function
(PGF); but no general solution is provided, and manual solution becomes
increasingly cumbersome as the number of priority levels grows beyond a small
number. No actual probability mass functions (PMFs) are computed. By contrast,
we present an explicit closed-form expression for the joint PGF given any
number of priority levels, that lends itself to practical numerical evaluation
of the joint PMF. In [12], explicit results are confined to the probability
that the system is empty, and the probability that a given queue is being
served at a random observation. In [13], the distribution of the maximum queue
length during a busy period in the presence of probabilistic feedback is also
given. The approach that we have adopted in the present paper may be viewed as
an extension of the method employed by Cohen [3] for the two-level case. With
it, we are able, in the words of Neuts [10], ‘to obtain actual results, that
is numbers and insight from numbers’.

Priority-level numbers significantly greater than two are encountered in
numerous real-world applications, such as health care [5]. In a hospital
emergency department (ED), arrivals are prioritized according to patient
acuity level [9]. Most hospitals operate with at least five acuity levels. In
a forthcoming paper, the results obtained here will be applied to the
ambulance ramping problem [1], in which arrivals to the ED by ambulance or as
walk-ins are categorized into three priority levels, corresponding to high,
intermediate or low patient acuity. The model is further complicated by the
fact that there are two arrival classes (ambulance and walk-in) each with
their own arrival rate and each of which contain patients of multiple priority
levels. Thus, there is an entanglement between arrival classes and priority
levels.

The number of servers is denoted by $c$ and the number of priority levels by
$K$. Each priority level is associated with a Poisson arrival rate
$\lambda_{\kappa}$, $\kappa=1,2,\ldots,K$, leading to a total arrival rate
$\lambda=\sum_{\kappa=1}^{K}\lambda_{\kappa}$. A common service rate $\mu$,
associated with an exponential distribution, is assumed for all priority
levels. Thus, the total traffic intensity is given by $r=\lambda/(c\mu)$. As
we are interested in the state-state queue-length distribution, possible
values of $r$ will be limited to the ergodic region $r<1$. The level traffic
intensity for priority $\kappa=1,2,\ldots,K$ is defined as
$r_{\kappa}=\lambda_{\kappa}/(c\mu)$, so that
$r=\sum_{\kappa=1}^{K}r_{\kappa}$. It is also convenient to introduce
priority-level fractions $0\leq\nu_{\kappa}\leq 1$, summing to unity, such
that $r_{\kappa}=r\nu_{\kappa}$. Alternatively, the $\nu_{\kappa}\geq 0$ may
be chosen without constraint, provided we set
$r_{\kappa}=r\nu_{\kappa}/\|\boldsymbol{\nu}\|_{1}$. This facilitates
exploration of different distributions among the priority levels for a given
constant total load on the system. To test the numerical performance of the
algorithms developed here, we shall fix $r$ and compute distributions for a
random sample of vectors $\boldsymbol{\nu}=(\nu_{1},\nu_{2},\ldots,\nu_{K})$
chosen independently from the unit probability simplex in $K$-dimensions
$\Delta_{K}\equiv\\{\boldsymbol{\nu}\in\mathbb{R}^{K}:\|\boldsymbol{\nu}\|_{1}=1,\boldsymbol{\nu}\geq
0\\}$. We also define the partial (per server) traffic intensity as
$\rho\equiv\lambda/\mu$ so that $r=c\rho$, in line with the notation of [8].

The rest of the paper is organized as follows: In Section 2, we set up the
stationary balance equations for the model. It is shown that, for any number
of priority levels, they can be stated in compact form comprising just a
single equation. This leads to a nearest-neighbour relationship among the
elements of the joint PMF, that is subsequently used as a diagnostic test of
the computational procedures developed later. It is also shown that the
balance equation can be solved directly for the joint PMF by means of a fixed-
point iteration (FPI). While the FPI does not constitute an efficient
computational algorithm, it provides an important benchmark for verifying the
correctness of the vastly more efficient methodology discussed in the
remainder of the paper. In Section 3, the general balance equation is used to
derive an explicit closed-form expression for the multi-variate PGF of the
joint distribution. It is shown that the correct marginals follow from the
joint PGF. Section 4 describes a multi-dimensional fast Fourier transform-
based (FFT) method that computes the joint PMF from the multi-variate PGF. It
is equally applicable to the marginal distributions. Measures of performance
for various diagnostic tests are presented in Section 5, the results of which
are discussed there. Conclusions follow in Section 6. Various technical
details appear in the Appendices.

##  2 Stationary Balance Equations

x Markov chain transitions for the two-level non-preemptive priority queue.

For ease of illustration, we shall begin by considering the concrete example
of a system with three priority levels: high, intermediate, low. Let the
single-subscripted quantity $p_{k}$ denote the probability that the queue is
empty and that $k$ of the $c$ servers are busy. Let the triple-subscripted
quantity $p_{\ell mn}$ denote the probability that there are $\ell$ high-
priority clients, $m$ intermediate-priority clients and $n$ low-priority
client in the queue, and that all servers are busy. Since, apart from the
identification $p_{c}=p_{000}$, the probabilities $p_{k}$ and $p_{\ell mn}$
are exhaustive and mutually exclusive, we must have that

$\sum_{k=0}^{c-1}p_{k}+\sum_{\ell,m,n=0}^{\infty}p_{\ell mn}=1\;.$ (1)

Figure 2 displays the state transition diagram for the two-level case, which
is easily drawn and captures all the qualitative features of the general
problem. With the definitions above, the stationary balance equations for the
three-level system are given by

$\displaystyle\begin{aligned} &(a)\quad&\lambda p_{k}&=(k+1)\mu
p_{k+1}&\text{for}\quad&0\leq k\leq c\;,\\\
&(b)\quad&(\lambda+c\mu)p_{000}&=\lambda
p_{c-1}+c\mu(p_{100}+p_{010}+p_{001})&&\\\
&(c)\quad&(\lambda+c\mu)p_{00n}&=\lambda_{3}p_{00n-1}+c\mu(p_{10n}+p_{01n}+p_{00n+1})&\text{for}\quad&n\geq
1\;,\\\
&(d)\quad&(\lambda+c\mu)p_{0mn}&=\lambda_{2}p_{0m-1n}+\lambda_{3}p_{0mn-1}+c\mu(p_{1mn}+p_{0m+1n})&\text{for}\quad&m\geq
1,n\geq 0\;,\\\ &(e)\quad&(\lambda+c\mu)p_{\ell
mn}&=\lambda_{1}p_{\ell-1mn}+\lambda_{2}p_{\ell m-1n}+\lambda_{3}p_{\ell
mn-1}+c\mu p_{\ell+1mn}&\text{for}\quad&\ell\geq 1,m,n\geq 0\;.\end{aligned}$
(2)

We adopt the convention that negative valued subscripts correspond to zero
probability. It follows immediately from (a) that

$p_{k}=\frac{1}{k!}\left(\frac{\lambda}{\mu}\right)^{k}p_{0}\;,$ (3)

for $k=0,1,\ldots,c$. Therefore,

$p_{000}\equiv p_{c}=\frac{1}{c!}\left(\frac{\lambda}{\mu}\right)^{c}p_{0}\;.$
(4)

With the observation that $\lambda p_{c-1}=c\mu p_{c}=c\mu p_{000}$, we can
recast the collection of balance equations above for the $p_{\ell mn}$ into
the single equation

$\begin{split}(\lambda+c\mu)p_{\ell
mn}&=\lambda_{1}p_{\ell-1mn}+\lambda_{2}p_{\ell m-1n}+\lambda_{3}p_{\ell
mn-1}\\\ &{}+c\mu\left(p_{\ell+1mn}+\delta_{0\ell}p_{\ell
m+1n}+\delta_{0\ell}\delta_{0m}p_{\ell
mn+1}+\delta_{0\ell}\delta_{0m}\delta_{0n}p_{\ell mn}\right)\;.\end{split}$
(5)

One should note that (5) is homogeneous in $p_{\ell mn}$, which means that its
solution is determined only up to an arbitrary scale factor. It is convenient
to set $\tilde{p}_{\ell mn}\equiv p_{\ell mn}/p_{000}$ so that
$\tilde{p}_{000}=1$. The normalization condition (1) then implies that

$\frac{1}{p_{000}}=\sum_{k=0}^{c-1}\frac{c!}{k!}\left(\frac{\mu}{\lambda}\right)^{c-k}+\sum_{\ell,m,n=0}^{\infty}\tilde{p}_{\ell
mn}\;.$ (6)

Thus, we may solve (5) for $\tilde{p}_{\ell mn}$ by setting
$\tilde{p}_{000}=1$, and subsequently determine $p_{000}$ from (6), so that
the probabilities can be recovered via $p_{\ell
mn}=p_{000}{\cdot}\tilde{p}_{\ell mn}$.

Next, let us write

$\Phi(\ell)\equiv\sum_{m,n}^{\infty}p_{\ell mn}\;,$ (7)

and observe that we also have

$\sum_{m,n=0}^{\infty}p_{\ell m-1n}=\sum_{m,n=0}^{\infty}p_{\ell
mn-1}=\Phi(\ell)\;.$ (8)

Moreover,

$\displaystyle\begin{aligned} \sum_{m,n=0}^{\infty}\delta_{0\ell}p_{\ell
m+1n}&=\delta_{0\ell}\Bigl{[}\Phi(0)-\sum_{n=0}^{\infty}p_{00n}\Bigr{]}\;,\\\
\sum_{m,n=0}^{\infty}\delta_{0\ell}\delta_{0m}p_{\ell
mn+1}&=\delta_{0\ell}\Bigl{[}\sum_{n=0}^{\infty}p_{00n}-p_{000}\Bigr{]}\;,\\\
\sum_{m,n=0}^{\infty}\delta_{0\ell}\delta_{0m}\delta_{0n}p_{\ell
mn}&=\delta_{0\ell}p_{000}\;,\end{aligned}$ (9)

and so we see that these three contributions sum to $\Phi(0)$. Accordingly,
(5) implies that

$(\lambda_{1}+c\mu)\Phi(\ell)=\lambda_{1}\Phi(\ell-1)+c\mu\left[\Phi(\ell+1)+\delta_{0\ell}\Phi(0)\right]\;,$
(10)

for $\ell=0,1,\ldots$. If we let
$\Delta\Phi(\ell)\equiv\Phi(\ell)-\Phi(\ell-1)$, then
$\Delta\Phi(\ell+1)=r_{1}\Delta\Phi(\ell)$ for $\ell=1,2,\ldots$, subject to
$\Delta\Phi(1)=-(1-r_{1})\Delta\Phi(0)$ and where $r_{1}=\lambda_{1}/(c\mu)$.
This is solved by $\Delta\Phi(\ell)=\varphi(\ell){\cdot}\Delta\Phi(0)$, with
$\varphi(\ell)\equiv-(1-r_{1})r_{1}^{\ell-1}$. It is equivalent to the first-
order recurrence equation

$\Phi(\ell)=\Phi(\ell-1)+\varphi(\ell){\cdot}\Phi(0)\;,$ (11)

for $\ell=1,2,\ldots$ which, in turn, is solved by

$\Phi(\ell)=\Phi(0){\cdot}\Bigl{[}1+\sum_{k=1}^{\ell}\varphi(k)\Bigr{]}=\Phi(0)r_{1}^{\ell}\;,$
(12)

for $\ell=0,1,\ldots$. The quantity $\Phi(\ell)$ represents the probability
that there are $\ell$ high-priority clients in the queue and that all servers
are busy. For $\ell\geq 1$, it simply gives the probability of $\ell$ high-
priority clients in the queue without qualification because having one or more
high-priority clients in the queue implies that all servers must be busy. Let
$\Psi(\ell)$ denote the probability of $\ell$ high-priority clients in the
queue. Then, for $\ell\geq 1$, $\Psi(\ell)=\Phi(\ell)=\Phi(0)r_{1}^{\ell}$ and

$1=\sum_{\ell=0}^{\infty}\Psi(\ell)=\Psi(0)+\sum_{\ell=1}^{\infty}\Phi(\ell)\;,$
(13)

from which it follows that $\Psi(0)=1-[r_{1}/(1-r_{1})]{\cdot}\Phi(0)$. Thus,
we may write

$\Psi(\ell)=\left[1-\Phi(0)/(1-r_{1})\right]\delta_{\ell
0}+\Phi(0)r_{1}^{\ell}\;.$ (14)

The normalization condition (1) implies that

$\sum_{\ell=0}^{\infty}\Phi(\ell)=1-P_{\text{NW}}\;,\quad
P_{\text{NW}}\equiv\sum_{k=0}^{c-1}p_{k}\;,$ (15)

where $P_{\text{NW}}$ represents the no-wait probability. We may equate this
with the summation over $\Phi(\ell)$ as given by (12) to obtain
$\Phi(0)=(1-r_{1})(1-P_{\text{NW}})$. Hence, (14) can be re-expressed in the
more direct form

$\Psi(\ell)=P_{\text{NW}}\delta_{\ell
0}+(1-P_{\text{NW}})(1-r_{1})r_{1}^{\ell}\;.$ (16)

This result applies equally to the general problem with an arbitrary number of
priority levels. The only undetermined component in the foregoing discussion
is the value of $p_{0}$, which represents the probability that the system is
empty, and is given by [8]

$\frac{1}{p_{0}}=\sum_{k=0}^{c}\frac{\rho^{k}}{k!}+\frac{\rho^{c+1}}{c!(c-\rho)}\;,$
(17)

where $\rho\equiv\lambda/\mu$. Knowledge of $p_{0}$ allows us to obtain
$p_{000}=p_{c}$. Another way of characterizing it is as follows: Let the
random variable $\mathcal{N}_{\text{sys}}$ represent the number of clients in
the system. It is well-known that, for some constant $A$,
$\Pr(\mathcal{N}_{\text{sys}}=n)=Ar^{n}$ for all $n\geq c$. By construction,
$p_{000}=\Pr(\mathcal{N}_{\text{sys}}=c)$. We also have
$P_{\text{NW}}=\Pr(\mathcal{N}_{\text{sys}}\leq c-1)$. Hence,

$\displaystyle\begin{aligned}
P_{\text{NW}}&=\sum_{n=0}^{c-1}\Pr(\mathcal{N}_{\text{sys}}=n)\\!\\!\\!&=1&-\sum_{n=c}^{\infty}\Pr(\mathcal{N}_{\text{sys}}=n)\\\
&=1-\frac{Ar^{c}}{1-r}&=1&-\frac{1}{1-r}\Pr(\mathcal{N}_{\text{sys}}=c)\;,\end{aligned}$
(18)

which leads to the result $p_{000}=(1-r)(1-P_{\text{NW}})$.

The stationary balance equations as given in (5) are easily extended to the
general problem of an arbitrary number $K$ of priority levels. In order to
establish a compact expression, we introduce the lattice vectors
$\mathbf{n}\equiv[n_{1},n_{2},\ldots,n_{K}]\in\mathbb{Z}^{K}$, and let
$\mathbf{e}_{\kappa}$, $\kappa=1,2,\ldots,K$ denote the standard unit
Cartesian coordinate basis vectors in $\mathbb{Z}^{K}$. Then the stationary
balance equations for $K$ priority levels can be expressed as

$(1+r)p_{\mathbf{n}}=\prod_{j=1}^{K}\delta_{0n_{j}}{\cdot}p_{\mathbf{n}}+\sum_{\kappa=1}^{K}\biggl{[}r_{\kappa}p_{\mathbf{n}-\mathbf{e}_{\kappa}}+\prod_{j=1}^{\kappa-1}\delta_{0n_{j}}{\cdot}p_{\mathbf{n}+\mathbf{e}_{\kappa}}\biggr{]}\;,$
(19)

where we recall the convention that $p_{\mathbf{n}}\equiv 0$ if $n_{\kappa}<0$
for any $\kappa=1,2,\ldots,K$. If we introduce the boundary set111Strictly
speaking, this is the boundary and beyond or, equivalently, the set of non-
interior points.

$\mathcal{B}\equiv\\{\mathbf{n}\in\mathbb{Z}^{K}:n_{\kappa}\leq 0\text{ for
some }\kappa=1,2,\ldots,K\\}\;,$ (20)

then, for all $\mathbf{n}\notin\mathcal{B}$, we have the interior stationary
balance equations

$p_{\mathbf{n}}=\frac{1}{1+r}\biggl{[}p_{\mathbf{n}+\mathbf{e}_{1}}+\sum_{\kappa=1}^{K}r_{\kappa}p_{\mathbf{n}-\mathbf{e}_{\kappa}}\biggr{]}\;.$
(21)

Thus every probability in the interior region is a positive weighted sum of
its lower nearest neighbours plus its upper highest priority neighbour. While
it is not a numerically stable proposition to attempt to first solve for
$p_{\mathbf{n}}$ on the boundary set $\mathcal{B}$ and then use (21) to
propagate the solution into the interior region, the relationship (21) does
provide a robust diagnostic test of where a candidate solution for
$p_{\mathbf{n}}$ behaves as it should. We shall call this the nearest-
neighbour test. The two-dimensional instance of this test was used extensively
for the two-level problem in [14]. The wait-conditional joint queue-length PMF
is given by $P(\mathbf{n})=(1-r)\tilde{p}_{\mathbf{n}}$. The full
unconditional joint PMF is then expressed as

$\displaystyle\begin{aligned}
P_{\text{full}}(\mathbf{n})&=P_{\text{NW}}\delta(\mathbf{n})+(1-P_{\text{NW}})P(\mathbf{n})\\\
&=P_{\text{NW}}\delta(\mathbf{n})+(1-P_{\text{NW}})(1-r)\tilde{p}_{\mathbf{n}}\;,\end{aligned}$
(22)

where $\delta(\mathbf{n})\equiv\prod_{\kappa=1}^{K}\delta_{0n_{\kappa}}$, and
recalling that $\tilde{p}_{\mathbf{0}}=1$.

The system (19) can be solved directly using an FPI. However, this
necessitates truncation of the problem to a maximum queue size for each
priority level. Consider the space $\mathcal{A}_{\infty}^{K}$ of countably
infinite matrices in $K$-dimensions, whose elements are non-negative and
indexed by the $K$-dimensional vector $\mathbf{n}=[n_{1},n_{2},\ldots,n_{K}]$,
with $n_{\kappa}=0,1,\ldots$, for each $\kappa=1,2\ldots,K$. In other words,
$\mathbf{n}\in{\mathbb{N}_{0}^{K}}$ — the space of $K$-tuples of non-negative
integers. We define the matrix mapping
$\mathcal{M}:\mathcal{A}_{\infty}^{K}\to\mathcal{A}_{\infty}^{K}$ by

$A^{\prime}=\mathcal{M}(A)\;,\quad
A^{\prime}_{\mathbf{n}}=\frac{1}{1+r}\left[\prod_{j=1}^{K}\delta_{0n_{j}}{\cdot}A_{\mathbf{n}}+\sum_{\kappa=1}^{K}\biggl{(}r_{\kappa}A_{\mathbf{n}-\mathbf{e}_{\kappa}}+\prod_{j=1}^{\kappa-1}\delta_{0n_{j}}{\cdot}A_{\mathbf{n}+\mathbf{e}_{\kappa}}\biggr{)}\right]\;.$
(23)

We observe that the mapping $\mathcal{M}$ is sum preserving:

$A^{\prime}=\mathcal{M}(A)\Rightarrow\sum_{\mathbf{n}\in{\mathbb{N}_{0}^{K}}}A^{\prime}_{\mathbf{n}}=\sum_{\mathbf{n}\in{\mathbb{N}_{0}^{K}}}A_{\mathbf{n}}\;,$
(24)

assuming that the sum is finite. This is equivalent to
$\|A^{\prime}\|_{1}=\|A\|_{1}$ for the $L_{1}$ matrix norm, and follows from
(10), which trivially extends to the general multi-level case.

We shall consider the collection of probabilities $p_{\mathbf{n}}$ as defining
a matrix $\mathsf{P}\in\mathcal{A}_{\infty}^{K}$. The $L_{1}$ and $L_{\infty}$
norms are given, respectively, by

$\|\mathsf{P}\|_{1}\equiv\sum_{\mathbf{n}}|\mathsf{P}_{\mathbf{n}}|\;,\quad\|\mathsf{P}\|_{\infty}\equiv\max_{\mathbf{n}}|\mathsf{P}_{\mathbf{n}}|\;.$
(25)

Both these norms exist and are finite for our matrix of (non-negative)
probabilities $\mathsf{P}$.

To solve the FPI directly, we must truncate the queue length in each dimension
(i.e. priority level) to a finite maximum size, $n_{\kappa}\leq
N_{\kappa}^{\text{max}}$ for each $\kappa=1,2,\ldots,K$. For the sake of
simplicity, we shall truncate in each of the $K$ dimensions to a common value
$N_{\text{max}}$. For the FPI, the matrix $\mathsf{P}$ will be initialized to
zero, except for $\mathsf{P}_{\mathbf{0}}=1$, on the grid
$[-1,N_{\text{max}}+1]^{K}$, but subsequently updated only on the sub-grid
$[0,N_{\text{max}}]^{K}$. Finite-size truncation will cause inevitable leakage
of probability at each step of the FPI. However, we can use the $L_{1}$
invariance to quantify the magnitude of the probability leakage, and adjust
for it by amortizing uniformly over all matrix elements. Also, after each FPI
step, we scale the matrix $\mathsf{P}$ to preserve
$\mathsf{P}_{\mathbf{0}}=1$. The correct normalization is applied at the end.
The details are provided in Algorithm 1.

Algorithm 1 Fixed-point iteration.

0: $(r_{1},\ldots,r_{K}),N_{\text{max}}$

0: $\mathsf{P}$

1: {Parameters:}

2: $\epsilon_{\text{tol}}=10^{-9}$

3: {Initialization:}

4: $\Delta=\infty$

5: $\mathsf{P}_{\mathbf{n}}\leftarrow\delta(\mathbf{n})$

6: while $\Delta>\epsilon_{\text{tol}}$ do

7: {Iteration:}

8: $\mathsf{P}^{\prime}=\mathcal{M}(\mathsf{P})$

9: {Probability leakage amortization:}

10: $p_{\text{leak}}=\|\mathsf{P}\|_{1}-\|\mathsf{P^{\prime}}\|_{1}>0$

11:
$\mathsf{P}^{\prime}\leftarrow\mathsf{P}^{\prime}+p_{\text{leak}}/(N_{\text{max}}+1)^{K}$

12: {Renormalization:}

13: $\mathsf{P}^{\prime}=\mathsf{P}^{\prime}/\mathsf{P}^{\prime}_{\mathbf{0}}$

14: {Convergence:}

15: $\Delta=\|\mathsf{P}^{\prime}-\mathsf{P}\|_{\infty}$

16: $\mathsf{P}\leftarrow\mathsf{P}^{\prime}$

17: end while

18: {Finalization:}

19: $\mathsf{P}\leftarrow(1-r){\cdot}\mathsf{P}$

Apart from the requirement that the queue length for each priority level be
truncated to a finite maximum value, the FPI method has other disadvantages,
as it suffers from the ‘curse of dimensionality’, and convergence becomes very
slow as the total traffic intensity $r$ gets close to unity. The slow
convergence in this region is exacerbated by the fact that the maximum queue
size needs to be set quite large, as probabilities of large queue sizes become
non-negligible. Its role in the present discussion is that it serves as a
benchmark for verifying the correctness for the much more efficient FFT method
developed in the following sections. Not only will it confirm that the joint
queue-length distribution has been computed correctly, there will be an
implied confirmation of the multi-variate PGF on which the FFT method rests.
Results of this comparison are illustrated in Figure (5.5) and discussed later
on in the results section. It is difficult to use Monte Carlo (MC) simulation
to check a multi-variate distribution in its entirety due to size requirements
and lack of appropriate hypothesis tests. We have, however, performed the more
straightforward task of checking marginal distributions against discrete-event
MC simulation. The FPI method does have the singular advantage that one may
expect it to be robust to the incorporation into the model of a variety of
complications (e.g. unequal service rates, or more complex queue disciplines),
whereas the closed-form PGF may not survive such modifications.

##  3 Probability Generating Function

Generalizing the approach of Cohen [3], let us introduce the collection of
multi-variate functions of $K-1$ continuous variables

$G_{\ell}(\mathbf{u})\equiv
P_{0}{\cdot}\sum_{m_{1},\ldots,m_{K-1}=0}^{\infty}\tilde{p}_{\ell m_{1}\cdots
m_{K-1}}u_{1}^{m_{1}}\cdots u_{K-1}^{m_{K-1}}\;,$ (26)

for $\ell=0,1,\ldots$, where $\mathbf{u}\equiv[u_{1},\ldots,u_{K-1}]$. With
the identification $P_{0}\equiv 1-r$, we have

$G_{0}(\mathbf{0})=1-r\;,\quad\sum_{\ell=0}^{\infty}G_{\ell}(\mathbf{1})=1\;,$
(27)

and $G_{\ell}(\mathbf{u})$ represents the PGF for the wait-conditional joint
queue-length PMF, such that

$P(\mathbf{n})=\prod_{\kappa=1}^{K-1}\frac{1}{n_{\kappa+1}!}\left.\frac{\partial^{n_{\kappa+1}}}{\partial
u^{n_{\kappa+1}}_{\kappa}}{\cdot}G_{n_{1}}(u_{1},\ldots,u_{K-1})\right|_{\mathbf{u}=\mathbf{0}}\;.$
(28)

Then, the PGF for the full unconditional PMF is given by

$G_{\text{full},\ell}(\mathbf{u})=P_{\text{NW}}\delta_{0\ell}+(1-P_{\text{NW}})G_{\ell}(\mathbf{u})\;.$
(29)

Summing over the stationary balance equations yields the linear recurrence
relations

$G_{\ell+1}+[\alpha-(1+r)]G_{\ell}+r_{1}G_{\ell-1}=0\;,$ (30)

for $\ell\geq 1$, where

$\alpha\equiv\mathbf{r}\cdot\mathbf{u}=\sum_{\kappa=1}^{K-1}r_{\kappa+1}u_{\kappa}\;.$
(31)

The characteristic equation reads

$\lambda^{2}+[\alpha-(1+r)]\lambda+r_{1}=0\;,$ (32)

and is solved by $\lambda=\lambda_{\pm}$ with

$\lambda_{\pm}=\tfrac{1}{2}\left[1+r-\alpha\pm\sqrt{(1+r-\alpha)^{2}-4r_{1}}\right]\;.$
(33)

When $\mathbf{u}=\mathbf{1}\equiv[1,\ldots,1]$, in which case
$\alpha=r-r_{1}$, we require, from the foregoing discussion of the high-
priority marginal, that $G_{\ell}\propto r_{1}^{\ell}$, and we have
$\lambda_{-}=r_{1}$ whereas $\lambda_{+}=1$. It follows immediately that

$G_{\ell}(\mathbf{u})=G_{0}(\mathbf{u})\lambda_{-}^{\ell}(\mathbf{u})\;,$ (34)

for $\ell=0,1,\ldots$. It remains to solve the problem for $\ell=0$.

For this purpose, let us write $\mathbf{n}\equiv(\ell,\mathbf{m})$ with
$\mathbf{m}=(m_{1},\ldots,m_{K-1})\in\mathbb{N}_{0}^{K-1}$ and
$q_{\mathbf{m}}^{\ell}\equiv p_{\ell\mathbf{m}}$. Then, for $\ell=0$, and on
setting

$Q_{\mathbf{m}}\equiv\sum_{j=1}^{K-1}q^{0}_{\mathbf{m}+\mathbf{e}_{j}}{\cdot}\prod_{i=1}^{j-1}\delta_{0m_{i}}\;,$
(35)

the birth-death equations (19) may be cast as

$(1+r)q^{0}_{\mathbf{m}}=\prod_{j=1}^{K-1}\delta_{j0m_{j}}{\cdot}q^{0}_{\mathbf{m}}+\sum_{j=1}^{K-1}r_{j+1}q^{0}_{\mathbf{m}-\mathbf{e}_{j}}+q^{1}_{\mathbf{m}}+Q_{\mathbf{m}}\;.$
(36)

Summing over the multi-indices $\mathbf{m}$ with the corresponding powers of
$\mathbf{u}$, we see that the stationary balance equations (19) require that

$(1+r)G_{0}(\mathbf{u})=G_{0}(\mathbf{0})+\alpha
G_{0}(\mathbf{u})+G_{1}(\mathbf{u})+\sum_{\mathbf{m}\in\mathbb{N}_{0}^{K-1}}Q_{\mathbf{m}}{\cdot}\prod_{j=1}^{K-1}u_{j}^{m_{j}}\;.$
(37)

Noting that $G_{1}(\mathbf{u})=\lambda_{-}G_{0}(\mathbf{u})$, and using the
identity $\lambda_{+}+\lambda_{-}=1+r-\alpha$, we arrive at

$\lambda_{+}G_{0}(\mathbf{u})=G_{0}(\mathbf{0})+\sum_{\mathbf{m}\in\mathbb{N}_{0}^{K-1}}Q_{\mathbf{m}}{\cdot}\prod_{j=1}^{K-1}u_{j}^{m_{j}}\;,$
(38)

which may be further manipulated to yield

$\displaystyle\begin{aligned} &\lambda_{+}G_{0}(u_{1},\ldots,u_{K-1})\\\
&=\sum_{k=1}^{K-1}\frac{1}{u_{k}}\left[G_{0}(0,\ldots,0,u_{k},\ldots,u_{K-1})-G_{0}(0,\ldots,0,u_{k+1},\ldots,u_{K-1})\right]+G_{0}(0,\ldots,0)\\\
&=\frac{1}{u_{1}}G_{0}(u_{1},\ldots,u_{K-1})+\sum_{k=2}^{K-1}\left(\frac{1}{u_{k}}-\frac{1}{u_{k-1}}\right)G_{0}(0,\ldots,0,u_{k},\ldots,u_{K-1})+\left(1-\frac{1}{u_{K-1}}\right)G_{0}(0,\ldots,0)\;.\end{aligned}$
(39)

On formally setting $u_{K}\equiv 1$, this may be simplified as

$G_{0}(u_{1},\ldots,u_{K-1})=\frac{1}{1/u_{1}-\lambda_{+}(u_{1},\ldots,u_{K-1})}\sum_{k=2}^{K}\left(\frac{1}{u_{k-1}}-\frac{1}{u_{k}}\right)G_{0}(0,\ldots,0,u_{k},\ldots,u_{K-1})\;.$
(40)

It should be noted that $G^{\prime}_{0}(u_{k},\ldots,u_{K-1})\equiv
G_{0}(0,\ldots,0,u_{k},\ldots,u_{K-1})$ is the result for an appropriately
aggregated $(K-k+1)$-dimensional problem, obtained by treating the $k$ highest
priority levels as a single high priority level with level traffic intensity
$r_{\text{hi}}=\sum_{\kappa=1}^{k}r_{\kappa}$. One may also note that
$G^{\prime}_{0}()=G_{0}(0,\ldots,0)=P_{0}$.

One immediate consequence of (40) is that the marginal distribution of the
aggregation of the top $p$ priority levels is geometric, for all
$p=1,2,\ldots,K$. From (34), the PGF for this distribution is

$G_{\text{agg}}^{(p)}(u)\equiv\sum_{\ell=0}^{\infty}u^{\ell}G_{\ell}(\mathbf{u}^{(p)})=\frac{G_{0}(\mathbf{u}^{(p)})}{1-u\lambda_{-}(\mathbf{u}^{(p)})}\;,$
(41)

where $\mathbf{u}^{(p)}\equiv(u{\cdot}\mathbf{1}_{p-1},\mathbf{1}_{K-p})$.
Only a single term survives the summation in (40) in the evaluation of
$G_{0}(\mathbf{u}^{(p)})$, to produce

$G_{\text{agg}}^{(p)}(u)=\frac{1-u}{[1-u\lambda_{+}(\mathbf{u}^{(p)})][1-u\lambda_{-}(\mathbf{u}^{(p)})]}G_{0}(\mathbf{0}_{p-1},\mathbf{1}_{K-p})\;.$
(42)

Now, we have

$\alpha^{(p)}\equiv\mathbf{r}\cdot\mathbf{u}^{(p)}=(\sigma_{p}-\sigma_{1})u+r-\sigma_{p}\;.$
(43)

where $\sigma_{p}\equiv\sum_{k=1}^{p}r_{k}$. It follows that we can write

$\lambda_{\pm}(\mathbf{u}^{(p)})=\tfrac{1}{2}\left[1+r_{\text{agg}}-r_{\text{lo}}u\pm\sqrt{(1+r_{\text{agg}}-r_{\text{lo}}u)^{2}-4r_{\text{hi}}}\right]\;,$
(44)

with

$r_{\text{lo}}=\sigma_{p}-\sigma_{1}\;,\quad r_{\text{hi}}=\sigma_{1}\;,\quad
r_{\text{agg}}=r_{\text{lo}}+r_{\text{hi}}=\sigma_{p}\;.$ (45)

Hence,

$\frac{1-u}{(1-u\lambda_{+})(1-u\lambda_{-})}=\frac{1}{1-\sigma_{p}u}\;,$ (46)

and the desired result follows, consistent with intuitive expectations.

The solution of the system (40) is derived in Appendix 7. Here, we explain the
overall strategy: We begin by noting that the recurrence relation (30) is
similar to equation (1.2.10) of Cohen [3], but with his analogue of
$r_{2}u_{1}$ replaced by our $\alpha$ defined in (31). The analogue of
(1.2.12) still holds, but the analogue of (1.2.11) acquires additional
complicating terms of the form $G_{0}(0,\ldots,0,u,v,w,\ldots)$. However,
these can be dealt with by means of a simple aggregation trick. This is
because the leading zeros mean that we looking at probabilities of various
numbers of low priority entities in the queue and no entities with the highest
$n$ priorities. But this is the same as no entities of a single highest
priority level aggregated from the highest $n$ priorities. Thus, we have to
solve a lower-dimensional problem. The relevant observation can be expressed,
for example in a 4-level problem, as $G_{0}(0,v,w)=G^{\prime}_{0}(v,w)$,
$G_{0}(0,0,w)=G^{\prime}_{0}(w)$, and of course
$G_{0}(0,0,0)=G^{\prime}_{0}()=P_{0}$. The PGFs with lesser numbers of
arguments refer to appropriately aggregated lower-dimensional problems.
Therefore, one has to start with Cohen’s two-dimensional problem and
successively work upwards to the full $K$-dimensional problem. The recurrence
relations take on a simpler, more transparent structure, that is amenable to
explicit solution, when one changes notation by writing the arguments of the
PGFs backwards. Accordingly, it is convenient to present the result for
$G_{0}(.)$ in terms of new variables $(z_{1},z_{2},\ldots,z_{K-1})$ where
ascending indices are associated with ascending priority level. Thus
$z_{k}=u_{K-k}$, in which case

$G_{0}(0,\ldots,0,u_{k},\ldots,u_{K-1})=G_{0}(0,\ldots,0,z_{K-k},\ldots,z_{2},z_{1})$
(47)

or, equivalently,
$G^{\prime}_{0}(u_{k},\ldots,u_{K-1})=G^{\prime}_{0}(z_{K-k},\ldots,z_{2},z_{1})$.
The result is then given for a new version of $G_{0}(.)$ re-parameterized
according to $G_{0}(z_{1},\ldots,z_{K-1})\equiv
G^{\prime}_{0}(z_{K-1},\ldots,z_{2},z_{1})$, which is just $G^{\prime}_{0}(.)$
with its arguments flipped.

A self-contained summary of Appendix 7 is as follows: To construct the PGF for
the problem of $K$ non-preemptive priority levels, we consider the joint PMF
$p_{\ell mn\cdots}$, where $\ell$ enumerates the number in the queue of the
highest priority level, and write the PGF as

$G_{\ell}(z_{1},\ldots,z_{K-1})\equiv\sum_{m_{1},\ldots,m_{K-1}=0}^{\infty}p_{\ell
m_{1}\cdots m_{K-1}}z_{K-1}^{m_{1}}\cdots z_{1}^{m_{K-1}}\;,$ (48)

so that $z_{1}$ is associated with the lowest priority level and $z_{K-1}$
with the next-to-highest priority level. We focus our attention on the
$\ell=0$ component $G_{0}(z_{1},\ldots,z_{K-1})$, and write $P_{0}\equiv
G_{0}(0,\ldots,0)$. Then, we have

$G_{0}(z_{1},\ldots,z_{K-1})=P_{0}{\cdot}\prod_{\kappa=1}^{K-1}\frac{1-z_{\kappa}\zeta_{+}(z_{1},\ldots,z_{\kappa-1})}{1-z_{\kappa}\zeta_{+}(z_{1},\ldots,z_{\kappa})}\;.$
(49)

The $\zeta_{\pm}(z_{1},\ldots,z_{\kappa})$ solve the quadratic equation

$\zeta^{2}+[\beta(z_{1},\ldots,z_{\kappa})-(1+r)]\zeta+\sigma_{K-\kappa}=0\;,$
(50)

for $\kappa=1,2,\ldots,K-1$. We have defined
$\sigma_{\kappa}\equiv\sum_{k=1}^{\kappa}r_{k}$, so that
$\sigma_{K-\kappa}=\sum_{k=1}^{K-\kappa}r_{k}$ and $\sigma_{K}=r$. Also,

$\beta(z_{1},\ldots,z_{\kappa})\equiv\sum_{k=1}^{\kappa}z_{k}r_{K+1-k}\;,$
(51)

for $\kappa=1,2,\ldots,K-1$. We adopt the conventions that $z_{0}\equiv 1$,
$\beta()\equiv 0$. The solutions $\zeta_{\pm}$ are given by

$\zeta_{\pm}(z_{1},\ldots,z_{\kappa})=\tfrac{1}{2}\left[1+r-\beta(z_{1},\ldots,z_{\kappa})\pm\sqrt{(1+r-\beta(z_{1},\ldots,z_{\kappa}))^{2}-4\sigma_{K-\kappa}}\right]\;.$
(52)

It follows that $\zeta_{+}()=1$, $\zeta_{-}()=r$.

The explicit expression for the two-level problem ($K=2$) is

$\displaystyle\begin{aligned}
G_{0}(z_{1})&=P_{0}{\cdot}\frac{1-z_{1}}{1-z_{1}\zeta_{+}(z_{1})}\\\
&=P_{0}{\cdot}\frac{1-z_{1}}{1-\dfrac{z_{1}}{2}\left[1+r-r_{2}z_{1}+\sqrt{(1+r-r_{2}z_{1})^{2}-4r_{1}}\right]}\;,\end{aligned}$
(53)

which is in agreement with Cohen [3], noting that $P_{0}=1-r$ if $G_{0}(.)$ is
to represent the wait-conditional PGF. The explicit expression for the three-
level problem ($K=3$) is

$\displaystyle\begin{aligned}
G_{0}(z_{1},z_{2})&=P_{0}{\cdot}\frac{1-z_{1}}{1-z_{1}\zeta_{+}(z_{1})}{\cdot}\frac{1-z_{2}\zeta_{+}(z_{1})}{1-z_{2}\zeta_{+}(z_{1},z_{2})}\\\
&=P_{0}{\cdot}\frac{1-z_{1}}{1-\dfrac{z_{1}}{2}\left[1+r-r_{3}z_{1}+\sqrt{(1+r-r_{3}z_{1})^{2}-4(r_{1}+r_{2})}\right]}\\\
&\quad{}\times\frac{1-\dfrac{z_{2}}{2}\left[1+r-r_{3}z_{1}+\sqrt{(1+r-r_{3}z_{1})^{2}-4(r_{1}+r_{2})}\right]}{1-\dfrac{z_{1}}{2}\left[1+r-r_{3}z_{1}-r_{2}z_{2}+\sqrt{(1+r-r_{3}z_{1}-r_{2}z_{2})^{2}-4r_{1}}\right]}\;.\end{aligned}$
(54)

The result for the full wait-conditional PGF is given by

$G_{\ell}(z_{1},\ldots,z_{K-1})=G_{0}(z_{1},\ldots,z_{K-1}){\cdot}\zeta_{-}^{\ell}(z_{1},\ldots,z_{K-1})\;,$
(55)

for $\ell=0,1,\ldots$.

###  3.1 Marginal Probabilities

We first consider the low-priority marginal PMF for the two-level ($K=2$)
problem. The low-priority marginal PGF is a function of a single continuous
complex variable $z$, and is obtained by summing (55) over the discrete high-
priority queue lengths, to yield

$G_{\text{lo}}(z)=\sum_{\ell=0}^{\infty}G_{\ell}(z)=\frac{(1-r)(1-z)}{[1-z\zeta_{+}(z)][1-\zeta_{-}(z)]}\;,$
(56)

with

$\zeta_{\pm}(z)=\tfrac{1}{2}\left[1+r-r_{2}z\pm\sqrt{(1+r-r_{2}z)^{2}-4r_{1}}\right]\;,$
(57)

which agrees with Cohen’s result [3]. By using the quadratic identities

$\zeta_{+}(z)+\zeta_{-}(z)=1+r-r_{2}z\;,\quad\zeta_{+}(z)\cdot\zeta_{-}(z)=r_{1}\;,$
(58)

which also imply that

$[1-\zeta_{+}(z)]{\cdot}[1-\zeta_{-}(z)]=r_{2}(z-1)\;,$ (59)

one is able to express $G_{\text{lo}}(z)$ in various equivalent and more
convenient forms:

$G_{\text{lo}}(z)=\frac{1-r}{r_{2}}{\cdot}\frac{r-\zeta_{-}(z)}{1-rz}=\frac{1-r}{\zeta_{+}(z)-r}\;.$
(60)

Since $\zeta_{+}(1/r)=r$, $\zeta_{-}(1/r)=r_{1}/r$, the first form shows that
$G_{\text{lo}}(z)$ has a pole at $z=1/r>1$. We also have $\zeta_{+}(1)=1$,
$\zeta_{-}(1)=r_{1}$, but it should be clear that the apparent pole at $z=1$
implied by (56) is spurious due to a cancellation with the numerator.

The $p$-th marginal for the general $K$-level problem, with
$p=1,2,\ldots,K-1$, has its PGF given by

$G_{\text{mrg}}^{(p)}(z_{p})\equiv\sum_{\ell=0}^{\infty}G_{\ell}(1,\ldots,z_{p},\ldots,1)=\frac{G_{0}(1,\ldots,z_{p},\ldots,1)}{1-\zeta_{-}(1,\ldots,z_{p},\ldots,1)}\;,$
(61)

and pertains to the $p$-th lowest priority level or, equivalently, the
$(K{+}1{-}p)$-th highest priority level. Each $G_{\text{mrg}}^{(p)}(z_{p})$
has the same functional form as $G_{\text{lo}}(z)$ for the two-level problem.
Let us relabel (57) as

$\zeta_{\pm}(z;r_{\text{hi}},r_{\text{lo}})=\tfrac{1}{2}\left[1+r_{\text{sum}}-r_{\text{lo}}z\pm\sqrt{(1+r_{\text{sum}}-r_{\text{lo}}z)^{2}-4r_{\text{hi}}}\right]\;,$
(62)

by setting $r_{\text{hi}}=r_{1}$, $r_{\text{lo}}=r_{2}$, $r_{\text{sum}}\equiv
r_{\text{hi}}+r_{\text{lo}}$, so that we also have

$G_{\text{lo}}(z;r_{\text{hi}},r_{\text{lo}})=\frac{1-r_{\text{sum}}}{\zeta_{+}(z;r_{\text{hi}},r_{\text{lo}})-r_{\text{sum}}}\;.$
(63)

Then
$G_{\text{mrg}}^{(p)}(z_{p})=G_{\text{lo}}(z;r_{\text{hi}},r_{\text{lo}})$ on
making the identifications

$r_{\text{lo}}=r_{K+1-p}\;,\quad
r_{\text{hi}}=\sum_{\kappa=1}^{K-p}r_{\kappa}\;.$ (64)

The derivation of this result is relegated to the second appendix. Its
interpretation is as follows: All priority levels above the $p$-th lowest
level are aggregated into a single high level whose level traffic intensity is
the sum of the aggregated values. The low level is identified with the $p$-th
lowest level and assigned its level traffic intensity, namely $r_{K+1-p}$.
Finally, all priority levels below the $p$-th are discarded, so the the total
traffic intensity for the equivalent two-level problem becomes
$r_{\text{sum}}=r_{\text{lo}}+r_{\text{hi}}$.

###  3.2 Exclusively-Low Probabilities

To consider the case where there are assumed to be only the lowest-priority
clients in the queue, we must set $z_{2}=z_{3}=\cdots=z_{K-1}=0$, to obtain
the PGF for the exclusively-low case

$G_{\text{xlo}}(z_{1})=G_{0}(z_{1},0,\ldots,0)\;,$ (65)

which is given by (49) and (55) as

$G_{\text{xlo}}(z_{1})=P_{0}{\cdot}\frac{1-z_{1}}{1-z_{1}\zeta_{+}(z_{1},0,\ldots,0)}\;.$
(66)

On setting $r_{\text{lo}}=r_{K}$,
$r_{\text{hi}}=\sigma_{K-1}=r-r_{\text{lo}}$, we can write

$\displaystyle\begin{aligned}
\beta(z_{1},\mathbf{0})&=r_{\text{lo}}z_{1}\;,\\\
\zeta_{\pm}(z_{1},\mathbf{0})&=\tfrac{1}{2}\left[1+r-r_{\text{lo}}z_{1}\pm\sqrt{(1+r-r_{\text{lo}}z_{1})^{2}-4r_{\text{hi}}}\right]\;.\end{aligned}$
(67)

Thus, we see that $\zeta_{\pm}(z_{1},\mathbf{0})$ is identical with its
counterpart for the two-level ($K=2$) problem, so that we can write
$\zeta_{\pm}(z_{1},\mathbf{0})=\zeta_{\pm}(z_{1})$. Using the identity

$\frac{1-z_{1}}{1-z_{1}\zeta_{+}(z_{1})}=\frac{1-\zeta_{-}(z_{1})}{\zeta_{+}(z_{1})-r}\;,$
(68)

we arrive at the result

$G_{\text{xlo}}(z_{1})=P_{0}{\cdot}\frac{1-\zeta_{-}(z_{1})}{\zeta_{+}(z_{1})-r}=P_{0}{\cdot}\left[1+r_{\text{lo}}{\cdot}\frac{z_{1}}{\zeta_{+}(z_{1})-r}\right]\;.$
(69)

This is equivalent to

$G_{\text{xlo}}(z_{1})=1-r+r_{\text{lo}}z_{1}G_{\text{lo}}(z_{1})\;.$ (70)

It follows that

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

where $P_{\text{lo}}(n)$ is the wait-conditional marginal PMF for the lowest-
priority level. One should note that this relationship is independent of the
number of priority levels $K$.

###  3.3 Exclusively-High Probabilities

To consider the case where there are assumed to be only the highest-priority
clients in the queue, we must set $z_{1}=z_{1}=\cdots=z_{K-1}=0$, to obtain
the PMF for the exclusively-high case

$P_{\text{xhi}}(\ell)=G_{\ell}(0,0,\ldots,0)=P_{0}\zeta_{-}^{\ell}(\mathbf{0})\;,$
(72)

as given by (49) and (55). On setting $r_{\text{hi}}=\sigma_{1}=r_{1}$,
$r_{\text{lo}}=r-r_{\text{hi}}$, we can write

$\displaystyle\begin{aligned} \beta(\mathbf{0})&=0\;,\\\
\zeta_{\pm}(\mathbf{0})&=\tfrac{1}{2}\left[1+r\pm\sqrt{(1+r)^{2}-4r_{\text{hi}}}\right]\;.\end{aligned}$
(73)

Therefore

$P_{\text{xhi}}(\ell)=(1-r){\cdot}\left[\frac{1+r-\sqrt{(1+r)^{2}-4r_{\text{hi}}}}{2}\right]^{\ell}\;,$
(74)

for $\ell=0,1,\ldots$. One should note that this result is independent of the
number of priority levels $K$.

##  4 FFT Mixture Method

x Error analysis for the FFT mixture method.

Let us begin by considering the numerical evaluation of the marginal PMFs.
Without loss of generality, we may restrict our attention to the low-priority
PMF for the two-level problem. We have

$P_{\text{lo}}(n)=\frac{1}{n!}\left.\frac{d^{n}}{dz^{n}}g(z)\right|_{z=0}\;,\quad
g(z)=\frac{1-r}{\zeta_{+}(z)-r}\;,$ (75)

where $g(z)=G_{\text{lo}}(z)$ as given by (60). Invoking Cauchy’s integral
theorem, we can write this as

$P_{\text{lo}}(n)=\oint_{\mathcal{C}}\frac{dz}{2\pi
i}\,\frac{g(z)}{z^{n+1}}\;,$ (76)

for $n=0,1,\ldots$, where the integration contour $\mathcal{C}$ encircles the
origin anti-clockwise and remains within the radius of convergence of the
Taylor expansion of $g(z)$, which is given by $\eta_{\text{c}}=1/r>1$. If we
choose $\mathcal{C}$ to be the circle of radius $\eta<\eta_{\text{c}}$ centred
about the origin, and write $z=\eta e^{-i\theta}$, for $0\leq\theta<2\pi$,
then

$P_{\text{lo}}(n)=\frac{1}{\eta^{n}}\int_{0}^{2\pi}\frac{d\theta}{2\pi}\
e^{in\theta}g\left(\eta e^{-i\theta}\right)\;.$ (77)

Approximation of the integral by an $N$-interval trapezoidal rule on the grid
$\theta_{k}=2\pi k/N$, $k=0,1,\ldots,N$, yields the finite sum

$P_{\text{lo}}(n)\simeq\frac{1}{\eta^{n}N}\sum_{k=0}^{N-1}e^{2\pi
ink/N}g\left(\eta e^{-2\pi ik/N}\right)\;,$ (78)

which is proportional to the inverse discrete Fourier transform of the
sequence $h(k)=g\left(\eta e^{-2\pi ik/N}\right)$, $k=0,1,\ldots,N-1$. When
$N$ is chosen to be a power of $2$, it can be implemented as an inverse fast
Fourier transform (IFFT), in which case we write
$P_{\text{lo}}(n)\simeq\eta^{-n}{\cdot}{\sf IFFT}[h](n)$. Extraction of queue-
length probabilities from a PGF by means of a discrete Fourier transform has
been previously considered by Daigle [4]. We shall adopt an alternative
general method due to Fornberg [7] that is better suited to the multi-variate
problem

The RHS of (77) is independent of the value of the contour radius $\eta$. This
is no longer true in (78) as a consequence of the approximation. The
dependence on the integration radius in using an FFT to numerically compute
derivatives of analytic functions has been studied by Bornemann [2]. He found
that there exists an optimal radius that minimizes the error, and that this
optimal value is close to the radius of convergence.

The dependence of the approximated solution on the contour radius can be used
to one’s advantage, as originally pointed out by Fornberg [7]. Consider the
numerical differentiation of the analytic function $g(z)$. According to the
FFT method, we must evaluate

$H(n;\eta)\equiv\eta^{-n}{\cdot}{\sf IFFT}_{k}[h(k;\eta)](n)\;,\quad
h(k,\eta)\equiv g\left(\eta e^{-2\pi ik/N}\right)\;,$ (79)

for some contour radius. We may equally evaluate the weighted mixture

$H_{M}(n)\equiv\sum_{m=1}^{M}f_{m}H(n;\eta_{m})\;,$ (80)

for some collection of contour radii $\eta_{m}$ and mixture coefficients
$f_{m}$ that sum to unity $\sum_{m=1}^{M}f_{m}=1$. As discussed in [7], the
first $M$ aliasing terms can be cancelled using (in principle) any $M$
distinct contour radii $\eta_{m}<\eta_{\text{c}}$ provided ones makes an
appropriate choice of the coefficients $f_{m}$. To achieve this, the mixture
coefficients $f_{m}$ must solve the matrix equation

$\begin{bmatrix}1&1&\cdots&1\\\
\eta_{1}^{N}&\eta_{2}^{N}&\cdots&\eta_{M}^{N}\\\
\eta_{1}^{2N}&\eta_{2}^{2N}&\cdots&\eta_{M}^{2N}\\\
\vdots&\vdots&\ddots&\vdots\\\
\eta_{1}^{(M-1)N}&\eta_{2}^{(M-1)N}&\cdots&\eta_{M}^{(M-1)N}\end{bmatrix}\\!\\!\begin{bmatrix}f_{1}\vphantom{1}\\\
f_{2}\vphantom{\eta_{M}}\\\ f_{3}\vphantom{\eta_{M}^{2}}\\\
\vdots\vphantom{\vdots}\\\
f_{M}\vphantom{\eta_{M}^{M-1}}\end{bmatrix}=\begin{bmatrix}1\\\ 0\\\ 0\\\
\vdots\\\ 0\end{bmatrix}\;.$ (81)

The solution for the mixture coefficients $f_{m}$ is given by

$\frac{1}{f_{m}}=\prod_{\begin{subarray}{c}\ell=1\\\ \ell\neq
m\end{subarray}}^{M}\left(1-\frac{\eta_{m}^{N}}{\eta_{\ell}^{N}}\right)\;,$
(82)

for $m=1,2,\dots,M$. Therefore, the residual errors after the FFT mixture are
given by

$\varepsilon_{M}(n)=c_{M\\!N}(n)\sum_{m=1}^{M}f_{m}\eta_{m}^{MN}=(-1)^{M-1}c_{M\\!N}(n)\prod_{m=1}^{M}\eta_{m}^{N}\;,$
(83)

for some constants $c_{M\\!N}(n)$. The result for the summation over $m$
follows directly from (82). One may note that, when applied to (75),
$c_{M\\!N}(n)=P_{\text{lo}}(n+M\\!N)$.

###  4.1 Marginal Distribution

Without loss of generality, we consider the marginal distribution for the low-
priority level in the two-level ($K=2$) problem. For any fixed $r$, worst case
performance of the FFT method is observed to occur as $r_{\text{hi}}\to
0^{+}$. Thus, we shall set $r_{\text{lo}}=r$, $r_{\text{hi}}=0$. In this case,
the wait-conditional PGF is trivially given by

$g(z)=\frac{1-r}{1-rz}=\sum_{n=0}^{\infty}p_{n}z^{n}\;,$ (84)

with $p_{n}=(1-r)r^{n}$, for $k=0,1,\ldots$. To approximate the $p_{n}$ via an
$N$-point FFT, we apply the trapezoidal rule on the grid $z_{k}=\eta e^{-2\pi
ik/N}$ for $k=0,1,\ldots,N-1$, to obtain

$p_{n}=\frac{1}{N}\sum_{k=0}^{N-1}g(z_{k})e^{2\pi ikn/N}\;.$ (85)

This leads us to consider the ratio

$\frac{p_{n}}{(1-r)r^{n}}=\frac{1}{N(r\eta)^{n}}\sum_{k=0}^{N-1}\frac{e^{2\pi
ikn/N}}{1-r\eta e^{-2\pi ik/N}}$ (86)

that will approach unity as $N\to\infty$. With $\xi\equiv r\eta<1$, we
introduce the ratio function

$R_{N}(\xi,n)\equiv\frac{1}{N\xi^{n}}\sum_{k=0}^{N-1}\frac{e^{2\pi
ikn/N}}{1-\xi e^{-2\pi ik/N}}\;.$ (87)

Given (86), we aim to achieve $R_{N}(\xi,n)\simeq 1$ for all
$n=0,1,\ldots,N-1$. We may also observe that

$R_{N}(\xi,n)=\frac{1}{\xi^{n}}I_{N}(\xi,n)=\frac{1}{\xi^{n}}\mathop{\sf
IFFT}\nolimits_{k}\left[(1-\xi e^{-2\pi ik/N})^{-1}\right](n)\;.$ (88)

The ratio function can be evaluated by explicitly summing the series to give

$R_{N}(\xi,n)=\frac{1}{1-\xi^{N}}\;,$ (89)

independent of $n$. Thus, $\delta R_{N}(\xi)\equiv
R_{N}(\xi,n)-1\simeq\xi^{N}$ for $\xi^{N}\ll 1$. If the desired discretization
relative error is set to $10^{-\alpha}$, then $\xi$ is determined by
$\xi^{N}=10^{-\alpha}$.

We shall be directly computing the quantity $I_{N}(\xi,n)$ in (88) via a
dedicated FFT routine. Since $R_{N}(\xi,n)\simeq 1$, we have that
$I_{N}(\xi,n)\sim\xi^{n}$. This cannot be too small for relevant values of $n$
due to machine arithmetic limits. In double-precision arithmetic, we must
ensure that $I_{N}(\xi,n)\gg 10^{-16}$ for the largest
$n=0,1,\ldots,N_{\text{max}}<N$ of interest, $N$ being the FFT size, which
must be greater than $N_{\text{max}}$ in order to prevent aliasing. If we take
the absolute FFT error (i.e. round-off error) to be
$\varepsilon_{\text{fft}}\simeq 10^{-15}$, then the relative FFT error is
$10^{-15}/\xi^{N_{\text{max}}}\sim 10^{-15+\alpha N_{\text{max}}/N}$. Best
numerical performance is obtained by equating the two competing and opposing
sources of error, i.e. the relative discretization error with the relative FFT
error. This yields the relationship

$\alpha=\frac{15}{1+N_{\text{max}}/N}$ (90)

for the overall relative error exponent. It implies that, for twelve decimal
places of accuracy ($\alpha=12$), we require that $N\simeq 4N_{\text{max}}$.

Let us now consider the FFT mixture scheme. In order to specify the $M$
contour radii, we introduce a spread parameter $s$, whose value we typically
take to be $s\simeq 0.05$, and spread factors $\varsigma_{m}\equiv
1-s{\cdot}(m-1)/(M-1)$, for $m=1,2,\ldots,M$. Then, we set
$\xi_{m}=\varsigma_{m}\xi$, which gives all of the radii in terms of the
largest one $\xi_{1}=\xi$, as yet undetermined. The smallest is given by
$\xi_{M}=\varsigma_{M}\xi=(1-s)\xi$. It is useful to observe that, in the
present context, (83) reads

$\varepsilon_{M}(n)/[(1-r)r^{n}]=(-1)^{M-1}\prod_{m=1}^{M}\xi_{m}^{N}\;,$ (91)

the RHS being independent of $n$.

We revisit the foregoing error analysis adapted to the mixture method, and
begin by estimating the FFT error. If the mixture ratio function as given by
the trapezoidal rule were the actual object of evaluation by means of an FFT,
then the exact result would be

$R^{\text{mix}}_{N_{\text{fft}}}(\boldsymbol{\xi},n)=\sum_{m=1}^{M}\frac{f_{m}(\boldsymbol{\xi})}{1-\xi_{m}^{N_{\text{fft}}}}\;,$
(92)

while the approximate result rendered by the FFT is

$\displaystyle\begin{aligned}
R^{\text{fft}}_{N_{\text{fft}}}(\boldsymbol{\xi},n)&=\sum_{m=1}^{M}\frac{f_{m}(\boldsymbol{\xi})}{1-\xi_{m}^{N_{\text{fft}}}}{\cdot}\left(1\pm\frac{\varepsilon_{\text{fft}}}{\xi_{m}^{n}}\right)\\\
&=R^{\text{mix}}_{N_{\text{fft}}}(\boldsymbol{\xi},n)\pm\varepsilon_{\text{fft}}\sum_{m=1}^{M}\frac{|f_{m}(\boldsymbol{\xi})|}{1-\xi_{m}^{N_{\text{fft}}}}{\cdot}\frac{1}{\xi_{m}^{n}}\;,\end{aligned}$
(93)

where $\varepsilon_{\text{fft}}$ denotes the absolute error in the FFT
computation due to machine arithmetic precision. Thus, the error in the ratio
function due to FFT precision is

$\displaystyle\begin{aligned} \Delta
R_{N_{\text{fft}}}^{\text{fft}}(\boldsymbol{\xi},n)&\equiv\left|R^{\text{fft}}_{N_{\text{fft}}}(\boldsymbol{\xi},n)-R^{\text{mix}}_{N_{\text{fft}}}(\boldsymbol{\xi},n)\right|\\\
&=\varepsilon_{\text{fft}}\sum_{m=1}^{M}\frac{f_{m}(\boldsymbol{\xi})}{1-\xi_{m}^{N_{\text{fft}}}}{\cdot}\frac{1}{\xi_{m}^{n}}\\\
&<\varepsilon_{\text{fft}}\sum_{m=1}^{M}\frac{|f_{m}(\boldsymbol{\xi})|}{1-\xi_{m}^{N_{\text{fft}}}}{\cdot}\frac{1}{\xi_{m}^{n}}\;,\end{aligned}$
(94)

where the final step represents a worst-case bound. Given that $\xi_{m}^{N}\ll
1$, we have

$\sum_{m=1}^{M}\frac{|f_{m}(\boldsymbol{\xi})|}{1-\xi_{m}^{N}}{\cdot}\frac{1}{\xi_{m}^{n}}\simeq\sum_{m=1}^{M}\frac{|f_{m}(\boldsymbol{\xi})|}{\xi_{m}^{n}}\simeq\frac{|f_{M}(\boldsymbol{\xi})|}{\xi_{M}^{n}}\;.$
(95)

So, we see that the FFT error is dominated by the smallest contour radius
$\xi_{M}$. It follows that it is estimated by

$\Delta
R_{N_{\text{fft}}}^{\text{fft}}(\boldsymbol{\xi},n)\sim\varepsilon_{\text{fft}}{\cdot}|f_{M}(\boldsymbol{\xi})|/\xi_{M}^{n}\sim\varepsilon_{\text{fft}}/\xi_{M}^{n}\;,$
(96)

since $|f_{M}(\boldsymbol{\xi})|\simeq 1$ for a sufficiently large spread $s$,
and we take $\varepsilon_{\text{fft}}=10^{-15}$.

To estimate the discretization error, we consider the mixture ratio function

$\displaystyle\begin{aligned}
R^{\text{mix}}_{N_{\text{fft}}}(\boldsymbol{\xi},n)&=\sum_{m=1}^{M}\frac{f_{m}(\boldsymbol{\xi})}{1-\xi_{m}^{N_{\text{fft}}}}\\\
&=\sum_{m=1}^{M}f_{m}(\boldsymbol{\xi})+\sum_{m=1}^{M}\xi_{m}^{N_{\text{fft}}}f_{m}(\boldsymbol{\xi})+\cdots\\\
&=1+(-1)^{M-1}(\xi_{1}\xi_{2}\cdots\xi_{M})^{N_{\text{fft}}}+\cdots\;.\end{aligned}$
(97)

Thus, the discretization error due to application of the trapezoidal rule is

$\Delta
R^{\text{tpz}}_{N_{\text{fft}}}(\boldsymbol{\xi},n)\equiv\left|R^{\text{mix}}_{N_{\text{fft}}}(\boldsymbol{\xi},n)-1\right|\simeq(g\xi)^{MN_{\text{fft}}}\;,$
(98)

where we have introduced the geometric mean

$g\equiv\biggl{(}\prod_{m=1}^{M}\varsigma_{m}\biggr{)}^{1/M}\;.$ (99)

We proceed to equate the trapezoidal-rule error of (98) with the FFT error of
(96) at the largest desired queue length $n=N_{\text{max}}$, and set them to a
common error level $10^{-\alpha}$, i.e.

$(g\xi)^{MN_{\text{fft}}}=10^{-15}/\xi_{M}^{N_{\text{max}}}=10^{-\alpha}\;.$
(100)

Then, after taking $\alpha=12$, $M=4$, eliminating $\xi$ from the equations
yields the relationship $N_{\text{fft}}=N_{\text{max}}/(1-N_{\text{max}}\chi)$
with $\chi\equiv\log_{10}((1-s)/g)$. For $s=0.05$, $M=4$, we obtain
$\chi=0.0037$. So, it suffices to choose $N_{\text{fft}}\simeq
N_{\text{max}}$. Consequently, the largest contour radius is determined as

$\xi=g^{-1}{\cdot}10^{-12/(MN_{\text{fft}})}\;,$ (101)

with $N_{\text{fft}}$ chosen to be the smallest power of two exceeding
$N_{\text{max}}$.

As a sanity check, we apply the FFT-mixture method as described here to the
computation of the low-priority marginal PMF for the two-level problem of the
case $r_{\text{hi}}=0$, and plot in Figure 4 various errors for the ratio
function as a function of queue length. The maximum queue length considered
was that sufficient to attain a tail probability $P_{\text{min}}=10^{-12}$.
The target error level was set to $10^{-12}$ (i.e. $\alpha=12$). The exact
(pre-discretization) ratio function is everywhere unity in the present case.
The overall error (blue curve) represents the difference in the computed ratio
function from unity. The FFT error (red curve) represents the difference in
the computed ratio function from the exact value of the discretized ratio
function. The discretization error (orange curve) represents the difference in
the exact discretized ratio function from unity. The results are consistent
with expectations.

###  4.2 Joint Distribution

The expression (49) for $G_{0}(\mathbf{z})$ contains spurious algebraic
singularities that cancel out between denominator and numerator, as can be
seen by observing that $\zeta_{+}(\mathbf{1})=1$. This has implications for
numerical evaluation, and limits the choice of integration contour radii in
the application of Cauchy’s theorem to $|z_{\kappa}|<1$ whereas, in principle,
one could extend this to at least $|z_{\kappa}|<1/r$. Fortunately, the
algebraic singularities can be eliminated by re-structuring the integrand into
an equivalent form according to the following argument: Let us write (49) as

$G_{0}(\mathbf{z})=P_{0}{\cdot}\prod_{\kappa=1}^{K-1}\mathscr{P}^{+}_{\kappa}(\mathbf{z})/\mathscr{Q}^{+}_{\kappa}(\mathbf{z})\;,$
(102)

where

$\mathscr{P}^{\pm}_{\kappa}(\mathbf{z})\equiv
1-z_{\kappa}\zeta^{\pm}_{\kappa-1}(\mathbf{z})\;,\quad\mathscr{Q}^{\pm}_{\kappa}(\mathbf{z})\equiv
1-z_{\kappa}\zeta^{\pm}_{\kappa}(\mathbf{z})\;,$ (103)

with
$\zeta^{\pm}_{\kappa}(\mathbf{z})\equiv\zeta_{\pm}(z_{1},\ldots,z_{\kappa})$,
for $\kappa=0,1,\ldots,K-1$. In particular, $\zeta^{+}_{0}()=1$,
$\zeta^{-}_{0}()=r$. Then, as derived in Appendix 9, we have the result

$\mathscr{P}^{+}_{\kappa}(\mathbf{z})/\mathscr{Q}^{+}_{\kappa}(\mathbf{z})=\mathscr{Q}^{-}_{\kappa}(\mathbf{z})/\mathscr{P}^{-}_{\kappa}(\mathbf{z})\;,$
(104)

which yields the equivalent form for the PGF

$G_{\ell}(z_{1},\ldots,z_{K-1})=P_{0}{\cdot}\prod_{\kappa=1}^{K-1}\frac{1-z_{\kappa}\zeta_{-}(z_{1},\ldots,z_{\kappa})}{1-z_{\kappa}\zeta_{-}(z_{1},\ldots,z_{\kappa-1})}{\cdot}\zeta_{-}^{\ell}(z_{1},\ldots,z_{K-1})\;.$
(105)

This result can also be established directly via an argument based on the fact
that the marginal queue-length distribution for the aggregation of the highest
$p$ priority levels is a geometric distribution, for all $p=1,2,\ldots,K$. The
details are presented in Appendix 10. It is also interesting to note that the
change of integration variables $\mathbf{z}\mapsto\mathbf{w}$ such that
$w_{\kappa}(z_{1},\ldots,z_{\kappa})=\zeta_{-}(z_{1},\ldots,z_{\kappa})$, for
$\kappa=1,2\ldots,K-1$, leads to the completely meromorphic form

$G_{\ell}(z_{1},\ldots,z_{K-1}){\cdot}\prod_{\kappa=1}^{K-1}dz_{\kappa}=P_{0}w^{\ell}_{K-1}{\cdot}\prod_{\kappa=1}^{K-1}\frac{1-z_{\kappa}w_{\kappa}}{1-z_{\kappa}w_{\kappa-1}}{\cdot}\prod_{\kappa=1}^{K-1}\frac{\sigma_{K-\kappa}-w^{2}_{\kappa}}{r_{\kappa+1}w^{2}_{\kappa}}dw_{\kappa}\;,$
(106)

where

$\displaystyle\begin{aligned}
z_{\kappa}&=\frac{w_{\kappa}(w^{2}_{\kappa-1}+\sigma_{K+1-\kappa})-\sigma_{K-\kappa}w_{\kappa-1}-w_{\kappa-1}w^{2}_{\kappa}}{r_{K+1-\kappa}w_{\kappa-1}w_{\kappa}}\\\
&=\frac{1}{r_{K+1-\kappa}}\left(w_{\kappa-1}+\frac{\sigma_{K+1-\kappa}}{w_{\kappa-1}}-w_{\kappa}-\frac{\sigma_{K-\kappa}}{w_{\kappa}}\right)\;,\end{aligned}$
(107)

and we adopt the convention that $z_{0}\equiv
r=\sum_{\kappa=1}^{K}r_{\kappa}$. The second product is the Jacobian of the
variable transformation. We have not yet found a practical use for this
representation.

Application of Cauchy’s integral theorem in each of $K{-}1$ dimensions yields
the expression for the joint PMF

$P(\ell,n_{K-1},\ldots,n_{1})=\prod_{\kappa=1}^{K-1}\biggl{\\{}\oint_{\mathcal{C}_{\kappa}}\frac{dz_{\kappa}}{2\pi
iz_{\kappa}^{n_{\kappa}+1}}\biggr{\\}}\,{\cdot}G_{\ell}(z_{1},\ldots,z_{K-1})\;,$
(108)

where the anti-clockwise closed contours around the origin
$\mathcal{C}_{\kappa}$ can be taken to be circles of radius
$\eta_{\kappa}<1/r$. For simplicity, we adopt a common radius
$\eta_{\kappa}=\eta$ in each dimension.

We introduce the multi-dimensional inverse FFT

$\mathscr{F}_{\ell}(\mathbf{n};\eta)=\eta^{-N(K-1)}{\cdot}{\sf
IFFT}_{\mathbf{k}}[G_{\ell}(\eta e^{-2\pi i\mathbf{k}/N})](\mathbf{n})\;,$
(109)

where, for simplicity, the FFT size $N$ is taken to be the same for every
dimension, and $\mathbf{n}=[n_{1},n_{2},\ldots,n_{K-1}]$. According to the
mixture method,

$P(\ell,n_{K-1},\ldots,n_{1})\simeq\sum_{m=1}^{M}f_{m}\mathscr{F}_{\ell}(\mathbf{n};\eta_{m})\;,$
(110)

with the contour radii $\eta_{m}$ and the coefficients $f_{m}$ chosen as for
the marginal distributions discussed in the preceding section.

Using the mixture, we execute $M$ runs of an FFT of size $N$ in each of the
$K{-}1$ dimensions. Without the mixture technique, comparable accuracy would
require a single FFT execution of size $MN$. Therefore the no-mixture/mixture
timing ratio for $K$ priority levels is

$\frac{(MN)^{K-1}\log_{2}(MN)^{K-1}}{M{\cdot}N^{K-1}\log_{2}(N^{K-1})}=M^{K-2}{\cdot}\left[1+\frac{\log_{2}M}{\log_{2}N}\right]\simeq
M^{K-2}\;.$ (111)

We see that there is no benefit in implementing the mixture method for the
two-level ($K=2$) problem, but it is increasingly advantageous as the number
of priority levels rises. Testing of numerous cases confirms this performance
differential.

##  5 Numerical Tests

###  5.1 Aggregation Test

As discussed for the two-level problem in [14], 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=\sum_{\kappa=1}^{K}r_{\kappa}$, which is known to be a simple geometric
distribution. Hence, the exact aggregate PMF is given by

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

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

$[\Delta\ln
P_{\text{agg}}](k)\equiv\ln(P_{\text{agg}}(k))-\ln(P_{\text{agg}}(k-1))=\ln
r\;,$ (113)

independent of $k$.

One diagnostic test of the 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 marginals
as only finite summations are required. The aggregate PMF is obtained from the
joint PMF as

$P_{\text{agg}}(k)=\sum_{\begin{subarray}{c}\mathbf{n}\in\mathbb{N}_{0}^{K}\\\
n_{1}+\cdots+n_{K}=k\end{subarray}}P(\mathbf{n})\;,$ (114)

for $k=0,1,2,\ldots$. We then consider the measure of performance (MOP)

$\Xi_{\text{agg}}\equiv-\max_{k\geq 1}\left\\{\log_{10}\left(|[\Delta\ln
P_{\text{agg}}](k)-\Delta\ln P_{\text{agg}}^{(\text{ex})}|\right)\right\\}\;,$
(115)

where the maximum is taken over all values $0\leq k\leq n_{\text{lim}}$ such
that $P_{\text{agg}}^{(\text{ex})}(k)>P_{\text{min}}>0$, for some threshold
level $P_{\text{min}}$ since one cannot expect the numerical methods to
maintain performance down to arbitrarily small tail probabilities. Since we
are working in double-precision arithmetic, 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.

###  5.2 Nearest-Neighbour Test

A direct consequence of the identity (21) is that the joint PMF at any given
interior point $\mathbf{n}$ is a positively weighted sum of the joint PMF
values at $K+1$ of its $2K$ nearest neighbours. We denote by
$P_{\text{nn}}(\mathbf{n})$ the PMF evaluated at the point $\mathbf{n}$ by
means of (21). Then, we consider the MOP

$\Xi_{\text{nn}}\equiv-\max_{\mathbf{n}\in\mathbb{N}^{K}}\left\\{\log_{10}\left(|\ln(P(\mathbf{n}))-\ln(P_{\text{nn}}(\mathbf{n}))|\right)\right\\}\;,$
(116)

where the maximum is taken over all values
$\mathbf{n}\in[1,n_{\text{lim}}]^{K}$ such that
$P(\mathbf{n})>P_{\text{min}}>0$.

###  5.3 Xhi-Test

We denote by $P_{\text{xhi}}^{(\text{ex})}(\ell)$ the theoretical exact
expression of the exclusively-high PMF as given by (74), which we compare with
the result obtained directly from the joint PMF as

$P_{\text{xhi}}(\ell)=P(\ell,0,\ldots,0)\;.$ (117)

If we define, for all $\ell>0$,

$[\Delta\ln
P_{\text{xhi}}](\ell)\equiv\ln(P_{\text{xhi}}(\ell))-\ln(P_{\text{xhi}}(\ell-1))\;,$
(118)

then $\Delta\ln
P_{\text{xhi}}^{(\text{ex})}=\ln\bigl([1+r-\sqrt{(1+r)^{2}-4r_{\text{hi}}}]/2\bigr{missing})$,
independent of $\ell$. The MOP for the xhi-test is then taken to be

$\Xi_{\text{xhi}}\equiv-\max_{\ell\geq 1}\left\\{\log_{10}\left(|[\Delta\ln
P_{\text{xhi}}](\ell)-\Delta\ln
P_{\text{xhi}}^{(\text{ex})}|\right)\right\\}\;,$ (119)

where the maximum is taken over all values $0\leq\ell\leq n_{\text{lim}}$ such
that $P_{\text{xhi}}^{(\text{ex})}(\ell)>P_{\text{min}}>0$. When applied to
the FFT method, the xhi-test simply checks how accurately the numerical FFT
mixture algorithm reproduces the Cauchy residue theorem for a simple pole at
the origin.

###  5.4 Xlo-Test

The exclusively-low PMF, as obtained directly from the computed joint PMF, is
given by

$P_{\text{xlo}}(n)=P(0,\ldots,0,n)\;.$ (120)

It is theoretically related to the lowest-priority marginal PMF according to
(71). This allows one to compare the computed $P_{\text{xlo}}(n)$ against a
known accurate evaluation of the marginal $P_{\text{lo}}(n)$. Such a marginal
PMF is provided by the quadratic recurrence method developed in [14] for the
two-level problem, recalling that the two-level problem is sufficient for the
calculation of all marginals of the multi-level problem. The relevant MOP is
taken to be

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

where the maximum is taken over all values $0<n\leq n_{\text{lim}}$ such that
$P_{\text{lo}}(n)>P_{\text{min}}>0$.

###  5.5 FPI Test

In the FPI test, we compute the joint PMF using both the FPI and FFT methods,
and check their agreement. The MOP is simply taken to be

$\Xi_{\text{fpi}}\equiv-\max_{\mathbf{n}\in\mathbb{N}^{K}}\left\\{\log_{10}\left(|\ln(P_{\text{fft}}(\mathbf{n}))-\ln(P_{\text{fpi}}(\mathbf{n}))|\right)\right\\}\;,$
(122)

where the maximum is taken over all values
$\mathbf{n}\in[1,n_{\text{lim}}]^{K}$ such that
$P_{\text{fft}}(\mathbf{n})>P_{\text{min}}>0$.

x Aggregation test for $K=4$ priority levels with $30$ randomized level
traffic intensities corresponding to total traffic intensity $r=0.9$. The
number of decimal places of agreement with the exact result is plotted on the
vertical axis as a function of aggregated queue length. The maximum aggregated
queue length of $100$ includes data points with PMF above
$P_{\text{min}}=2.4\times 10^{-6}$.

x Xlo-test for $K=4$ priority levels with $30$ randomized level traffic
intensities corresponding to total traffic intensity $r=0.9$. The number of
decimal places of agreement is plotted on the vertical axis as a function of
lowest-priority queue length. The maximum queue length of $60$ includes data
points with PMF above $P_{\text{min}}=1.0\times 10^{-6}$.

x Xhi-test for $K=4$ priority levels with $30$ randomized level traffic
intensities corresponding to total traffic intensity $r=0.9$. The number of
decimal places of agreement is plotted on the vertical axis as a function of
highest-priority queue length. The maximum queue length of $100$ includes data
points with PMF above $P_{\text{min}}=1.0\times 10^{-20}$.

x Nearest-neighbour test for $K=4$ priority levels with $30$ randomized level
traffic intensities corresponding to total traffic intensity $r=0.9$. The
number of decimal places of agreement is plotted on the vertical axis as a
function of lowest-priority queue length. All joint-PMF data points above
$P_{\text{min}}=1.0\times 10^{-10}$ that occur within a maximum queue length
of $100$ in each priority dimension were considered.

x FFT versus FPI test for $K=3$ priority levels with $30$ randomized level
traffic intensities corresponding to total traffic intensity $r=0.75$. The
number of decimal places of agreement is plotted on the vertical axis as a
function of lowest-priority queue length. All joint-PMF data points above
$P_{\text{min}}=1.0\times 10^{-10}$ that occur within a maximum queue length
of $100$ in each priority dimension were considered.

###  5.6 Results

Figures 5.5–5.5 present the results of the numerical tests described above.
MOP values prior to worst-case minimization, relevant to the FFT mixture
computations, are displayed on the vertical axes against the relevant queue
lengths. Each test pertains to $K=4$ priority levels repeated with $30$
randomized level traffic intensities, all corresponding to total traffic
intensity $r=0.9$, as shown by the blue lines and dots. The red curve picks
out an individual run. The cyan curve delineates the overall worst-case
behaviour. The PMF array is generated up to a maximum queue length of
$N_{\text{max}}=100$ in each priority dimension. We see that
$\Xi_{\alpha}>9.5$ for $\alpha=\text{agg},\text{xhi},\text{xlo},\text{nn}$.

Figure 5.5 compares the joint PMF arrays as computed via the FFT and FPI
methods, for $K=3$ priority levels with $30$ randomized level traffic
intensities summing to $r=0.75$. All elements whose probabilities exceeded a
threshold tail level of $P_{\text{min}}>10^{-10}$ were included. The worst
case for each lowest-priority queue length is plotted. We see that
$\Xi_{\text{fpi}}\gtrsim 2$. Decreasing performance the further one progresses
into the tail is a reflection of the truncation error in the FPI method.

Finally, in Figure 5.6, we plot the results of the combined set of tests for
$K=7$ priority levels with a randomly generated set of level traffic
intensities that sum to a total traffic intensity of $r=0.9$. The PMF array is
generated up to a maximum queue length of $N_{\text{max}}=15$ in each priority
dimension. For the nearest-neighbour test, array elements with probabilities
above $P_{\text{min}}=10^{-6}$ are considered, and the results are presented
as the worst case for each lowest-priority queue length. The other curves are
self-explanatory. It can be observed that the algorithm holds up well in this
large example, where the FFT size is 0.25 GB, and the size of the PMF array is
2 GB.

The results obtained confirm (i) numerical stability, (ii) internal
consistency, (iii) consistency between methods, and (iv) agreement with
theoretical exact values.

x All FFT tests for $K=7$ priority levels with a random set of level traffic
intensities corresponding to total traffic intensity $r=0.9$. The number of
decimal places of agreement is plotted on the vertical axis as a function of
the relevant queue length. Joint-PMF data points were computed to a maximum
queue length of $15$ in each priority dimension. For the nearest-neighbour
test, all all points above $P_{\text{min}}=1.0\times 10^{-6}$ were considered.

##  6 Conclusions

Simple methods for accurate computation of the joint queue-length distribution
for a non-preemptive multi-level priority queue have been developed. An
explicit multi-variate PGF for the joint queue-length distribution has been
derived in closed form for the first time. With the aid of an FFT mixture
method, it is used to numerically compute the joint queue-length PGF and the
queue-length marginals. A direct method based on a fixed-point iteration of
the stationary balance equations has also been developed. Comparable explicit
and exact results are not available elsewhere in the existing literature.

##  7 Multi-Variate PGF

In order to solve (40) for

$G^{\prime}_{0}(u_{k},\ldots,u_{K-1})=G_{0}(0,\ldots,0,u_{k},\ldots,u_{K-1})\;,$
(123)

we consider a general problem for a family of functions $G_{0}(.)$ that are
differentiated from one another according to the number of arguments that they
take. The system to be solved has the following structure: For each
$n=1,2,\ldots,K-1$,

$G_{0}(x_{1},\ldots,x_{n})=\frac{1}{1/x_{1}-\lambda_{+}(x_{1},\ldots,x_{n})}\sum_{k=2}^{n+1}\left(\frac{1}{x_{k-1}}-\frac{1}{x_{k}}\right)G_{0}(x_{k},\ldots,x_{n})\;,$
(124)

with $x_{n+1}\equiv 1$, and $G_{0}()=P_{0}$. The functions
$\lambda_{\pm}(x_{1},\ldots,x_{n})$ are solutions of the quadratic equation

$\zeta^{2}-[1+r-\alpha(x_{1},\ldots,x_{n})]\zeta+\sigma_{K-n}=0\;,$ (125)

with

$\alpha(x_{1},\ldots,x_{n})=\sum_{k=1}^{n}x_{k}r_{K-n+k}\;.$ (126)

We observe that
$G^{\prime}_{0}(u_{k},\ldots,u_{K-1})=G_{0}(x_{1},\ldots,x_{n})$, being the
application with $n=K-k$ arguments, in which case

$\alpha(u_{k},\ldots,u_{K-1})=\sum_{\ell=1}^{K-k}u_{\ell+k-1}r_{\ell+k}=\sum_{\ell=k}^{K-1}u_{\ell}r_{\ell+1}\;,$
(127)

which involves only the $K-k$ lowest priorities with level traffic intensities
$r_{k+1},\ldots,r_{K}$, while the constant term of the quadratic
$\sigma_{k}=\sum_{\ell=1}^{k}r_{\ell}$ aggregates the disjoint highest $k$
priorities with level traffic intensities $r_{1},\ldots,r_{k}$. This is
consistent with the interpretation of $G^{\prime}_{0}(u_{k},\ldots,u_{K-1})$
as solving the $(K{+}1{-}k)$-level sub-problem with the $k$ highest priorities
aggregated into a single level while keeping the total traffic $r$ intensity
constant.

Next, we introduce

$\displaystyle\begin{aligned} J(x_{n},\ldots,x_{1})&\equiv
G_{0}(x_{1},\ldots,x_{n})\;,\\\
\zeta_{\pm}(x_{n},\ldots,x_{1})&\equiv\lambda_{\pm}(x_{1},\ldots,x_{n})\;,\end{aligned}$
(128)

and set $z_{k}\equiv x_{n+1-k}$ for $k=1,2,\ldots,n$, so that $z_{1}=x_{n}$,
$z_{n}=x_{1}$, and the convention $x_{n+1}=1$ translates to $z_{0}=1$. Then
(124) becomes

$\displaystyle\begin{aligned}
J(z_{1},\ldots,z_{n})&=\frac{1}{1/z_{n}-\zeta_{+}(z_{1},\ldots,z_{n})}\sum_{k=2}^{n+1}\left(\frac{1}{z_{n+2-k}}-\frac{1}{z_{n+1-k}}\right)J(z_{1},\ldots,z_{n+1-k})\\\
&=\frac{1}{1/z_{n}-\zeta_{+}(z_{1},\ldots,z_{n})}\sum_{\ell=0}^{n-1}\left(\frac{1}{z_{\ell+1}}-\frac{1}{z_{\ell}}\right)J(z_{1},\ldots,z_{\ell})\;,\end{aligned}$
(129)

under the change of summation variable $\ell=n+1-k$, and where
$J()=G_{0}(0,\ldots,0)=P_{0}$. Let us now write $J_{\ell}(\mathbf{z})\equiv
J(z_{1},\ldots,z_{\ell})$ and $J_{0}\equiv J()$, so that the subscript $\ell$
indicates the number of elements (i.e. dimension) of the vector argument
$\mathbf{z}$. Then, we have

$J_{n}=\frac{1}{1/z_{n}-\zeta_{+}(z_{1},\ldots,z_{n})}\sum_{\ell=0}^{n-1}\left(\frac{1}{z_{\ell+1}}-\frac{1}{z_{\ell}}\right)J_{\ell}\;,$
(130)

for $n=1,2,\ldots,K-1$. Finally, let
$J^{\prime}_{n}\equiv\left(1/z_{n}-\zeta_{+}(z_{1},\ldots,z_{n})\right)J_{n}$
so that we obtain the recurrence

$J^{\prime}_{n}=\sum_{\ell=0}^{n-1}W_{\ell}J^{\prime}_{\ell}\;,$ (131)

with

$W_{\ell}=\left(\frac{1}{z_{\ell+1}}-\frac{1}{z_{\ell}}\right)\frac{1}{1/z_{\ell}-\zeta_{+}(z_{1},\ldots,z_{\ell})}\;,$
(132)

for $\ell\geq 1$, and $W_{0}=1/z_{1}-1$.

For $n\geq 1$, we can write

$J^{\prime}_{n+1}=W_{n}J^{\prime}_{n}+\sum_{\ell=0}^{n-1}W_{\ell}J^{\prime}_{\ell}=(1+W_{n})J^{\prime}_{n}\;,$
(133)

which is solved by

$J^{\prime}_{n}=\prod_{\ell=1}^{n-1}(1+W_{\ell}){\cdot}J^{\prime}_{1}\;,\quad
J^{\prime}_{1}=W_{0}J^{\prime}_{0}\;,$ (134)

where we have set $J^{\prime}_{0}=P_{0}$. After some algebraic manipulation,
this leads to the explicit representation

$J_{n}(\mathbf{z})=P_{0}{\cdot}\prod_{\ell=1}^{n}\frac{1-z_{\ell}\zeta_{+}(z_{1},\ldots,z_{\ell-1})}{1-z_{\ell}\zeta_{+}(z_{1},\ldots,z_{\ell})}\;,$
(135)

where the identity $\zeta_{+}()=1$ has been invoked. The quantity that is
ultimately of interest for the $K$-level problem is
$J_{K-1}(\mathbf{z})=G_{0}(\mathbf{z})=G_{0}(z_{1},\ldots,z_{K-1})$, where the
arguments of $G_{0}(.)$ are now reversed.

To complete the specification of the solution, we recall that

$\alpha(u_{m},\ldots,u_{K-1})=\sum_{k=m}^{K-1}u_{k}r_{k+1}\;,$ (136)

and set $z_{m}\equiv u_{K-m}$ for $m=1,2\ldots,K-1$, in order to define the
function $\beta(\mathbf{z})$ via

$\beta(u_{K-1},\ldots,u_{m})\equiv\alpha(u_{m},\ldots,u_{K-1})\;,$ (137)

so that

$\beta(z_{1},\ldots,z_{K-m})=\sum_{k=m}^{K-1}z_{K-k}r_{k+1}\;.$ (138)

Equivalently, since $m$ is arbitrary,

$\beta(z_{1},\ldots,z_{n})=\sum_{\ell=1}^{n}z_{\ell}r_{K+1-\ell}\;.$ (139)

The functions $\zeta_{\pm}(z_{1},\ldots,z_{n})$ solve the quadratic equation

$\zeta^{2}-[1+r-\beta(z_{1},\ldots,z_{n})]\zeta+\sigma_{K-n}=0\;.$ (140)

The two branches are given by

$\zeta_{\pm}(z_{1},\ldots,z_{n})=\tfrac{1}{2}\left[1+r-\beta(z_{1},\ldots,z_{n})\pm\sqrt{(1+r-\beta(z_{1},\ldots,z_{n}))^{2}-4\sigma_{K-n}}\right]\;,$
(141)

for $n=1,2,\ldots,K-1$. Noting that $\sigma_{K}=\sum_{k=1}^{K}r_{k}=r$, and
with the convention $\beta()\equiv 0$, we recover the relationships
$\zeta_{+}()=1$, $\zeta_{-}()=r$ from extension to the case $n=0$.

##  8 Extraction of Marginals

In this appendix, we show that all marginals for the general $K$-level
priority problem can be inferred from knowledge of the structure of the PGF
for the two-level problem. It is convenient to introduce the notation

$\displaystyle\begin{aligned}
\xi^{\pm}_{k}(\mathbf{z})&\equiv\zeta_{\pm}(z_{1},\ldots,z_{k})\;,\\\
\beta_{k}(\mathbf{z})&\equiv\beta(z_{1},\ldots,z_{k})\;,\end{aligned}$ (142)

where the integer subscript indicates the number of elements (i.e. the
dimension) of the vector argument. Then, from (61) and (105), it follows that

$\sum_{\ell=0}^{\infty}G_{\ell}(\mathbf{z})=\sum_{k=0}^{K-1}\frac{1-z_{k}\zeta^{-}_{k}(\mathbf{z})}{1-z_{k+1}\zeta^{-}_{k}(\mathbf{z})}\;,$
(143)

where we have set $z_{K}\equiv 1$, and noted that

$P_{0}=1-r=1-\zeta_{-}()=1-\zeta^{-}_{0}(\mathbf{z})\;.$ (144)

To extract the $p$-th marginal, we must consider the vector argument
$\mathbf{z}=\mathbf{z}^{(p)}$ where we set $z_{k}=1$ for all $k\neq p$. Thus
$\mathbf{z}^{(p)}\equiv\mathbf{1}+(z_{p}-1)\mathbf{e}_{p}$. We may observe
that the $k$-th term of the product on the RHS of (143) is equal to unity,
expect when $k=p-1$ or $k=p$. This leads to a representation of the PGF for
the $p$-th marginal as

$G_{\text{mrg}}^{(p)}(z_{p})=\sum_{\ell=0}^{\infty}G_{\ell}(\mathbf{z}^{(p)})=\frac{1-\zeta^{-}_{p-1}(\mathbf{z}^{(p)})}{1-z_{p}\zeta^{-}_{p-1}(\mathbf{z}^{(p)})}{\cdot}\frac{1-z_{p}\zeta^{-}_{p}(\mathbf{z}^{(p)})}{1-\zeta^{-}_{p}(\mathbf{z}^{(p)})}\;,$
(145)

for $p=1,2,\ldots,K-1$. Since the identity
$\beta_{k}(\mathbf{1})=r-\sigma_{K-k}$ implies that
$\zeta^{+}_{k}(\mathbf{1})=1$, $\zeta^{-}_{k}(\mathbf{1})=\sigma_{K-k}$, we
have

$\zeta^{-}_{p-1}(\mathbf{z}^{(p)})=\zeta^{-}_{p-1}(\mathbf{1})=\sigma_{K+1-p}\;.$
(146)

Next, we observe that

$\beta_{p}(\mathbf{z}^{(p)})=\beta_{p}(1,\ldots,z_{p},\ldots,1)=\sum_{k=1}^{p-1}r_{K+1-k}+r_{K+1-p}z_{p}\;.$
(147)

Setting

$r_{\text{lo}}=r_{K+1-p}\;,\quad
r_{\text{hi}}=\sum_{k=1}^{K-p}r_{k}=\sigma_{K-p}\;,\quad r_{\text{sum}}\equiv
r_{\text{lo}}+r_{\text{hi}}=\sigma_{K+1-p}\;,$ (148)

we obtain

$r-\beta_{p}(\mathbf{z}^{(p)})=r_{\text{sum}}-r_{\text{lo}}z_{p}\;.$ (149)

Therefore,

$\displaystyle\begin{aligned}
\zeta^{\pm}_{p}(\mathbf{z}^{(p)})&=\tfrac{1}{2}\left[1+r-\beta_{p}(\mathbf{z}^{(p)})\pm\sqrt{\left(1+r-\beta_{p}(\mathbf{z}^{(p)})\right)^{2}-4\sigma_{K-p}}\right]\\\
&=\tfrac{1}{2}\left[1+r_{\text{sum}}-r_{\text{lo}}z_{p}\pm\sqrt{(1+r_{\text{sum}}-r_{\text{lo}}z_{p})^{2}-4r_{\text{hi}}}\right]\\\
&=\zeta_{\pm}(z_{p})\;,\end{aligned}$ (150)

where the final identification with $\zeta_{\pm}(z)$ for the two-level problem
is associated with the mappings $r_{1}\leftarrow r_{\text{hi}}$,
$r_{2}\leftarrow r_{\text{lo}}$, $r\leftarrow r_{\text{sum}}$. Consequently,

$G_{\text{mrg}}^{(p)}(z_{p})=\frac{1-r_{\text{sum}}}{1-r_{\text{sum}}z_{p}}{\cdot}\frac{1-z_{p}\zeta_{-}(z_{p})}{1-\zeta_{-}(z_{p})}\;,$
(151)

for $p=1,2,\ldots,K-1$. Hence, on comparing with (105) for the case $K=2$, we
see that the desired result for the marginal PGFs is established.
Alternatively, by multiplying numerator and denominator by
$1-\zeta_{+}(z_{p})$ and using the identities given in (58), one recovers the
first form in (60).

##  9 PGF Ratio

Considering the quantities $\mathscr{P}^{\pm}_{\kappa}(\mathbf{z})$,
$\mathscr{Q}^{\pm}_{\kappa}(\mathbf{z})$ as defined in (103), we shall show
that

$\mathscr{P}^{+}_{\kappa}(\mathbf{z})\mathscr{P}^{-}_{\kappa}(\mathbf{z})=\mathscr{Q}^{+}_{\kappa}(\mathbf{z})\mathscr{Q}^{-}_{\kappa}(\mathbf{z})\;,$
(152)

which implies that one can eliminate the removable singularities in (102) by
making the substitution

$\mathscr{P}^{+}_{\kappa}(\mathbf{z})/\mathscr{Q}^{+}_{\kappa}(\mathbf{z})\mapsto\mathscr{Q}^{-}_{\kappa}(\mathbf{z})/\mathscr{P}^{-}_{\kappa}(\mathbf{z})\;.$
(153)

In order to first evaluate $\mathscr{P}^{+}_{\kappa}\mathscr{P}^{-}_{\kappa}$,
we let $a\equiv 4\sigma_{K+1-\kappa}$ and $b\equiv 1+r-\beta_{\kappa-1}$. Then
it is easy to see that

$\displaystyle\begin{aligned}
\mathscr{P}^{+}_{\kappa}\mathscr{P}^{-}_{\kappa}&=1-bz_{\kappa}+\frac{a}{4}z^{2}_{\kappa}\\\
&=1-(1+r-\beta_{\kappa-1})z_{\kappa}+\sigma_{K+1-\kappa}z^{2}_{\kappa}\;.\end{aligned}$
(154)

Similarly, if we now let $a^{\prime}\equiv 4\sigma_{K-\kappa}$ and
$b^{\prime}\equiv 1+r-\beta_{\kappa}$, then

$\displaystyle\begin{aligned}
\mathscr{Q}^{+}_{\kappa}\mathscr{Q}^{-}_{\kappa}&=1-b^{\prime}z_{\kappa}+\frac{a^{\prime}}{4}z^{2}_{\kappa}\\\
&=1-(1+r-\beta_{\kappa})z_{\kappa}+\sigma_{K-\kappa}z^{2}_{\kappa}\\\
&=1-(1+r-\beta_{\kappa-1}-r_{K+1-\kappa}z_{\kappa})z_{\kappa}+\sigma_{K-\kappa}z^{2}_{\kappa}\\\
&=1-(1+r-\beta_{\kappa-1})z_{\kappa}+(\sigma_{K-\kappa}+r_{K+1-\kappa})z^{2}_{\kappa}\\\
&=1-(1+r-\beta_{\kappa-1})z_{\kappa}+\sigma_{K+1-\kappa}z^{2}_{\kappa}\;.\end{aligned}$
(155)

The desired result follows from the equality of the right-hand sides of the
foregoing pair of equations.

##  10 Alternative Derivation of the Joint PGF

Recalling from (55) that
$G_{\ell}(z_{1},\ldots,z_{K-1})=G_{0}(z_{1},\ldots,z_{K-1})\zeta_{-}^{\ell}$,
we obtain for the full $K$-dimensional joint PGF,

$G(z_{1},\ldots,z_{K})\equiv\sum_{\ell=0}^{\infty}z_{K}^{\ell}G_{\ell}(z_{1},\ldots,z_{K-1})=\frac{G(z_{1},\ldots,z_{K-1},0)}{1-z_{K}\zeta_{K-1}(\mathbf{z})}\;,$
(156)

where $\zeta_{K-1}(\mathbf{z})\equiv\zeta_{-}(z_{1},\ldots,z_{K-1})$ such that
the subscript denotes the number of arguments. We introduce the notation
$G(z_{1},\ldots,z_{p}|z)\equiv G(z_{1},\ldots,z_{p},z,\ldots,z)$ for any
$z\in\mathbb{C}$ and $p=0,1,\ldots,K$. For $p=0$, we have
$G({\cdot}|z)=G(z,\ldots,z)$, while for $p=K$, we have
$G(z_{1},\ldots,z_{K}|z)=G(z_{1},\ldots,z_{K})$, independent of $z$.

Since the marginal distribution of the aggregation of the top $p$ priority
levels is also geometric, as argued in the discussion below (41), we have that

$G(z_{1},\ldots,z_{K-p}|z)=\frac{G(z_{1}\ldots,z_{K-p}|0)}{1-z\zeta_{K-p}(\mathbf{z})}\;,$
(157)

for arbitrary $z$ and $p=1,2,\ldots,K$, where the rates
$\zeta_{K-p}(\mathbf{z})$ are to be determined. When $p=1$, we recover (156).
Setting $z=z_{K-p+1}$ in (157) yields

$G(z_{1},\ldots,z_{K-p}|z_{K-p+1})=\frac{G(z_{1}\ldots,z_{K-p}|0)}{1-z_{K-p+1}\zeta_{K-p}(\mathbf{z})}\;,$
(158)

for $p=1,\ldots,K$. Next we write (157) as

$G(z_{1},\ldots,z_{K-p+1}|z)=\frac{G(z_{1}\ldots,z_{K-p+1}|0)}{1-z\zeta_{K-p+1}(\mathbf{z})}\;,$
(159)

and choose $z_{K-p+1}=z$ so that

$G(z_{1},\ldots,z_{K-p}|z)=\frac{G(z_{1}\ldots,z_{K-p+1}|0)}{1-z\zeta_{K-p+1}(z_{1},\ldots,z_{K-p},z)}\;.$
(160)

Now we set $z=z_{K-p+1}$ to obtain

$G(z_{1},\ldots,z_{K-p}|z_{K-p+1})=\frac{G(z_{1}\ldots,z_{K-p+1}|0)}{1-z_{K-p+1}\zeta_{K-p+1}(\mathbf{z})}\;.$
(161)

We observe that (158) and (161) express equations for the same quantity. Thus,
on equating the RHS of each equation, we arrive at the recursion

$G(z_{1},\ldots,z_{k}|0)=\frac{1-z_{k}\zeta_{k}(\mathbf{z})}{1-z_{k}\zeta_{k-1}(\mathbf{z})}{\cdot}G(z_{1},\ldots,z_{k-1}|0)\;,$
(162)

for $k=1,2,\ldots,K-1$, which leads to the result

$G_{0}(z_{1},\ldots,z_{K-1})=G(z_{1},\ldots,z_{K-1}|0)=P_{0}\prod_{k=1}^{K-1}\frac{1-z_{k}\zeta_{k}(\mathbf{z})}{1-z_{k}\zeta_{k-1}(\mathbf{z})}\;,$
(163)

where we have used the fact that $P_{0}=G({\cdot}|0)=G(0,\ldots,0)$.

In order to determine the rates $\zeta_{k}(\mathbf{z})$, we multiply the
stationary balance equation (19) by $z_{1}^{n_{K}}\cdots
z_{K-k}^{n_{k+1}}{\cdot}z^{n_{k}+\cdots+n_{1}}$, and sum over all indices
$n_{1},\ldots n_{K}$ to obtain

$\left[\beta_{k}(\mathbf{z})-r+\frac{\sigma_{K-k}}{z\zeta_{k}(\mathbf{z})}\right]G(z_{1},\ldots,z_{k}|z)=G(z_{1},\ldots,z_{k}|0)\;,$
(164)

where

$\beta_{k}(\mathbf{z})=\sum_{\kappa=1}^{k}z_{\kappa}r_{K-k+1}\;,\quad\sigma_{k}=\sum_{\kappa=1}^{k}r_{\kappa}\;,$
(165)

and we have used the relation

$\left.\frac{d}{dz}G(z_{1},\ldots,z_{k}|z)\right|_{z=0}=\zeta_{k}(\mathbf{z})G(z_{1},\ldots,z_{k}|0)\;,$
(166)

that follows from (157). We also use (157), expressed as

$G(z_{1},\ldots,z_{k}|0)=(1-z\zeta_{k}(\mathbf{z}))G(z_{1},\ldots,z_{k}|z)\;,$
(167)

to obtain

$\left[z^{2}\zeta^{2}_{k}-(1+r-\beta_{k})z\zeta_{k}+\sigma_{K-k}\right]G(z_{1},\ldots,z_{k}|z)=0\;,$
(168)

provided $z\neq 0$. Setting $z=1$ yields

$\zeta^{2}_{k}-(1+r-\beta_{k})\zeta_{k}+\sigma_{K-k}=0\;,$ (169)

which is solved by (52). Therefore, the result (105) is established.

## Acknowledgments

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

## References

  * Almehdawe et al. [2013] Almehdawe E, Jewkes B, Q-M He (2013) A Markovian queueing model for ambulance offload delays. _European Journal of Operational Research_ 226(3):602–614.
  * Bornemann [2011] Bornemann F (2011) Accuracy and stability of computing higher-order derivatives of analytic function by cauchy integrals. _Foundations of Computational Mathematics_ 11(1):1–63.
  * 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.
  * Daigle [1989] Daigle J (1989) Queue length distributions from probability generating functions via discrete Fourier transforms. _Operations Research Letters_ 8(4):229–236.
  * Elalouf and Wachtel [2022] Elalouf A, Wachtel G (2022) Queueing problems in emergency departments: A review of practical approaches and research methodologoies. _Operations Research Forum_ 3(1):2:1–2:46.
  * Elmelegy [2010] Elmelegy A (2010) _An Exact Solution for the Problem of M/M/ $C$/$K$ Non-Preemptive Priority Queue Using State Equilibrium Equations_. Ph.D. thesis, Illinois Institute of Technology, Chicago, IL, USA, UMI Number: 3435818\.
  * Fornberg [1981] Fornberg B (1981) Numerical differentiation of analytic functions. _ACM Transactions on Mathematical Software_ 7(4):542–547.
  * Gnedenko and Kovalenko [1989] Gnedenko B, Kovalenko I (1989) _Introduction to Queueing Theory_ (Boston, MA, USA: Birkhäuser), second edition.
  * 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.
  * Neuts [1984] Neuts M (1984) Matrix-analytic methods in queuing theory. _European Journal of Operational Research_ 15(1):2–12.
  * 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.
  * Wignall [1973] Wignall T (1973) Priority queuing system with and without feedback. _Operations Research_ 21(3):764–776.
  * Wignall and Enns [1973] Wignall T, Enns E (1973) The joint stationary multivariate queue length distribution in a single server queueing system with n queues, arbitrary priorities, and a general probabilistic inter-queue transition matrix. _Management Science_ 19(7):778–782.
  * Zuk and Kirszenblat [2023] Zuk J, Kirszenblat D (2023) Explicit results for the distributions of queue lengths for a non-preemptive two-level priority queue, URL http://dx.doi.org/10.2139/ssrn.4574550, submitted for publication to _Annals of Operations Research_.


# A Generalized Weighted Loss for SVC and MLP

Filippo Portera

###### Abstract

Usually standard algorithms employ a loss where each error is the mere
absolute difference between the true value and the extrapolation, in case of a
regression task. In the present, we introduce several error weighting schemes
that are a generalization of the consolidated routine. We study both a binary
classification model for Support Vector Classification and a regression net
for Multy-layer Perceptron. Results proves that the error is never worse than
the standard procedure and several times it is better.

Machine Learning, Binary Classification, SVC, Regression, MLP.

##  1 Introduction

We would like to show that a standard loss generalization for binary
classification (in our case we have chosen SVC and MLP) could produce results
not worse w.r.t. the consolidated loss. In fact, the possibility that a given
data-set presents non-IID samples can be exploited by these generalized
losses.

The loss studied to generalize SVC and the full optimization problem are:

$P=\frac{1}{2}||\vec{v}||^{2}+C\sum_{i=1}^{l}\xi_{i}w_{i}$ (1)

suject to:

$y_{i}(\vec{v}^{\prime}\vec{x_{i}}+b)\geq 1-\xi_{i}\hskip
28.45274pti\in{1,..,l}$ (2)

and:

$\xi_{i}\geq 0\hskip 28.45274pti\in{1,..,l}$ (3)

where $\vec{v}$ represents the linear weights of the extraploator function,
$l$ is the number of training examples, $C$ is a trade-off hyper-parameter,
$\xi_{i}$ is the error on sample $i$, and $w_{i}$ are some scalar weights that
are a function of a distribution $s_{i}$ of the samples:

$s_{i}=\sum_{j=1}^{l}(e^{-\gamma_{S}||\vec{x}_{i}-\vec{x}_{j}||^{2}})$ (4)

Other distribution can be adopted (e.g., 1 + the RBF norm instead of the RBF
dot product). And let:

$sy_{i}=y_{i}y_{j}\sum_{j=1}^{l}(e^{-\gamma_{S}||\vec{x}_{i}-\vec{x}_{j}||^{2}})$
(5)

with $\gamma_{S}$ additional hyper-parameter. Here lies the complexity of the
algorithm since this calculation is $O(l^{2})$. Perhaps it can be overtaken
with pattern sampling or, in the case of MLP with a sort of weights learning.

$w_{i}=f(s_{i})$ (6)

This implies a quadratic problem that is different from traditional SVC:

The Lagrangian would be:

$L=\frac{1}{2}||\vec{v}||^{2}+C\sum_{i=1}^{l}\xi_{i}w_{i}-\sum_{i=1}^{l}\alpha_{i}(y_{i}(\vec{v}^{\prime}\vec{x_{i}}+b)-1+\xi_{i})$
(7)

$-\sum_{i=1}^{l}\eta_{i}\xi_{i}$ (8)

subject to:

$\alpha_{i}\geq 0\hskip 28.45274pti\in{1,..,l}$ (9)

$\eta_{i}\geq 0\hskip 28.45274pti\in{1,..,l}$ (10)

Applying the KKT condition for optimaility:

$\frac{\partial
L}{\partial\vec{v}}=\vec{v}-\sum_{i=1}^{1}\alpha_{i}y_{i}\vec{x_{i}}=0\Rightarrow\vec{v}=\sum_{i=1}^{1}\alpha_{i}y_{i}\vec{x_{i}}$
(11)

$\frac{\partial
L}{\partial\vec{\xi}}=\vec{\alpha}+\vec{\eta}-C\vec{w}=0\Rightarrow\vec{\alpha}\leq
C\vec{w}$ (12)

$\frac{\partial L}{\partial b}=\sum_{i=1}^{l}\alpha_{i}y_{i}=0$ (13)

Thus, the dual becomes:

$D=\sum_{i=1}^{l}\alpha_{i}-\frac{1}{2}\sum_{i=1}^{l}\sum_{j=1}^{l}\alpha_{i}y_{i}\alpha_{j}y_{j}K(\vec{x_{i}},\vec{x_{j}})$
(14)

subject to:

$0\leq\alpha_{i}\leq Cw_{i}\hskip 28.45274pti\in{1,..,l}$ (15)

$\sum_{i=1}^{l}\alpha_{i}y_{i}=0$ (16)

This is very similar to standard SVC dual [4], apart the constraints on the
lagrangian multipliers.

We wrote an ad-hoc quadratic optimizer for this problem111The code of this
work is available at OSF GWL Project, with a SMO-like method ([2]).

We iteratively select 2 distinct multipliers and we modify them with an
attempt to improve the dual objective function:

$\alpha_{i}^{t+1}=\alpha_{i}^{t}+\nu y_{i}$ (17)

$\alpha_{j}^{t+1}=\alpha_{j}^{t}-\nu y_{j}$ (18)

The motivation is the enforcement of the second dual constraint on the
$\sum_{i=1}^{l}\alpha_{i}y_{i}=0$.

The $\nu$ in the optimal direction is obtained deriving $D$ by $\nu$ as it has
been done in section 5.1 of ([5]).

This direction is:

$\nu=\frac{y_{j}-y_{i}-\sum_{p=1}^{l}\alpha_{p}y_{p}K(\vec{x_{j}},\vec{x_{p}})+\sum_{p=1}^{l}\alpha_{p}y_{p}K(\vec{x_{i}},\vec{x_{p}})}{K(\vec{x_{i}},\vec{x_{i}})-2K(\vec{x_{i}},\vec{x_{j}})+K(\vec{x_{j}},\vec{x_{j}})}$
(19)

Once the candidate $\nu$ has been determined, it has to be clipped in order to
satisfy the constraints on both the multipliers.

At each iteration we compute $b$ with the suport vectors that lie in the
margin (for which $0<\alpha_{i}<Cw_{i}$) as it has been reported in How to
calculate $b$.

The kernel used to compute $K(\vec{x},\vec{y})$ is RBF with hyper-parameter
$\gamma_{K}$. The whole procedure is iterated $50l^{2}$ times for each
training problem.

##  2 Related works

In ([1]) they learn the loss weights directly from the training and validation
sets. They assert that there is a substantial improvement in the
generalization error and they also provide theoretical bounds.

##  3 Method

We use the acronym GWL for Generalized Weighted Loss.

We tried 4 distinct algorithms: the Python 3 package sklearn.svm.SVC, GWL SVC
with $w_{i}=1$, GWL (here we mean the generalized loss with $w_{i}$’s built as
described), and GWL with random weights. We would like to know if, in the
general case, the optimal solutions use $w_{i}$ not equal to $1$. We have
selected at least 8 cases of study, to determine the weight $w_{i}$ of a
sample $i$. Therefore, some evaluated weighting functions are:

  1. 1.

$w_{i}=\sqrt{s_{i}}$ (20)

  2. 2.

$w_{i}=s_{i}$ (21)

  3. 3.

$w_{i}=s_{i}^{2}$ (22)

  4. 4.

$w_{i}=\frac{1}{\sqrt{s_{i}}}$ (23)

  5. 5.

$w_{i}=\frac{1}{s_{i}}$ (24)

  6. 6.

$w_{i}=\frac{1}{s_{i}^{2}}$ (25)

  7. 7.

$w_{i}=sy_{i}$ (26)

  8. 8.

$w_{i}=1+\textrm{rand}[0,1]$ (27)

The case $8$ is useful to show that a weighting scheme based on the training
distribution is more convenient w.r.t. a random weighting scheme.

##  4 Results

We explored a 2 dimensional hyper-parameters grid for sklearn.svm.svc,
involving $\gamma_{K}$ and $C$. While we used the additional hyper-parameter
$\gamma_{S}$ to generete the loss weights. That is the reason why the
experiments with loss weights take more time to terminate. Obviously, the
second grid is an extension of (it covers) the first one. Those are the
results for the 5-fold cross-validation with data-sets extracted from the UCI
website, and opportunely treated (double or inconsistent samples removed,
shuffling):

Algorithm | Data-set | Mean F1 | Time
---|---|---|---
sklearn.svm.SVC | Ionosphere | 0.968638 | 0m3,130s
GWL SVC | Ionosphere | 0.970651 | 62m34,796s
GWL(1) | Ionosphere | 0.977172 | 175m44,060s
GWL(2) | Ionosphere | 0.977172 | 175m41,022s
GWL(3) | Ionosphere | 0.977172 | 187m1,032s
GWL(4) | Ionosphere | 0.977172 | 187m59,323s
GWL(5) | Ionosphere | 0.977359 | 3h:07m:41s
GWL(6) | Ionosphere | 0.977538 | 2h:58m:42s
GWL(8) | Ionosphere | 0.977292 | 3h:24m:32s
GWL(8) | Ionosphere | 0.977292 | 3h:08m:31s
GWL(8) | Ionosphere | 0.974767 | 3h:07m:36s
GWL(8) | Ionosphere | 0.975011 | 3h:05m:40
Table 1: SVC and GWL with Ionosphere data-set Algorithm | Data-set | Mean F1 | Time
---|---|---|---
sklearn.svm.SVC | Sonar | 0.886610 | 0m2,396s
GWL SVC | Sonar | 0.904489 | 16m19,019s
GWL(1) | Sonar | 0.909337 | 28m8,391s
GWL(2) | Sonar | 0.916513 | 31m58,852s
GWL(3) | Sonar | 0.908717 | 36m26,726s
GWL(4) | Sonar | 0.913580 | 40m51,567s
GWL(5) | Sonar | 0.916303 | 45m46,03s
GWL(6) | Sonar | 0.916671 | 41m:17,92s
GWL(7) | Sonar | 0.911098 |
GWL(8) | Sonar | 0.907057 |
Table 2: SVC and GWL with Sonar data-set Algorithm | Data-set | Mean F1 | Time
---|---|---|---
sklearn.svm.SVC | Breast | 0.959825 | 0m3,387s
GWL SVC | Breast | 0.958628 | 174m47,138s
GWL(1) | Breast | 0.963625 | 432m3,673s
GWL(2) | Breast | 0.963896 | 448m14,994s
GWL(3) | Breast | 0.967909 | 443m59,346s
GWL(4) | Breast | 0.966109 | 401m18,642s
GWL(5) | Breast | 0.964666 | 6h:48m:01s
GWL(6) | Breast | 0.964666 | 6h:34m:13s
GWL(8) | Breast | 0.961837 | 8h:02m:57s
Table 3: SVC and GWL with Breast data-set Algorithm | Data-set | Mean F1 | Time
---|---|---|---
sklearn.svm.SVC | Statlog | 0.610351 | 0m38,297s
GWL SVC | Statlog | 0.644108 | 1725m29,269s
GWL(1) | Statlog | 0.651278 | 5370m49,360s
GWL(2) | Statlog | 0.651278 | 5499m18,497s
GWL(3) | Statlog | 0.644329 | 5638m2,544s
GWL(4) | Statlog | 0.649529 | 5692m55,211s
GWL(5) | Statlog | 0.651947 | 5351m40,975s
GWL(6) | Statlog | 0.652786 | 5546m26,466s
GWL(8) | Statlog | 0.646312 | 83h:50m:08s
Table 4: SVC and GWL with Statlog data-set

We also have tried 2 MLP nets with PyTorch on a regression task with
$w_{i}=s_{i}$ and results are interesting (but the random initialization of
the net weights should be considered in this case: Wine data-set: 3961
samples, 11 features; MLP: 100, 50, 20, 1, and Wine data-set, different MLP
100, 80, 40, 1 nodes per layer. The theory underneath deep neural
architectures can be foun in [3].

Standard MLP MAE | gamma Best | Best Loss MLP MAE
---|---|---
0.5871212 | 10 |  0.57575756 (1)
0.5694444 | 10 |  0.5530303 (2)
0.5568182 | 1 |  0.53661615 (3)
0.5580808 | 0.01 |  0.540404 (4)
0.54924244 | 0.1 |  0.54671717 (5)
0.510101 | 100 | 0.510101 (6)
Table 5: First MLP with wine data-set Standard MLP MAE | gamma Best | Best Loss MLP MAE
---|---|---
0.5694444 | 10 |  0.5580808 (1)
0.53661615 | 1 |  0.5290404 (2)
0.5378788 | 100 | 0.5378788 (3)
0.5770202 | 10 |  0.5580808 (4)
0.54924244 | 0.01 |  0.5252525 (5)
0.56565654 | 1 |  0.5555556 (6)
Table 6: Second MLP with wine data-set

In this scenario it would be useful to determine the difference between eah
couple of values, to understand which is the strategy that, in the most of the
cases, performs best with the test set. An idea is to learn weights, starting
to run in parallel $n$ nets with different random weight vectors and selecting
at each parallel one the best vector in terms of MAE and perturbing it and re-
run the procedure for a given amount of iterations.

GWL has been written in C. The regression code for the wine data-set has been
written in Python 3.10 and torch.

Hardware employed: a notebook with 8 cores Intel(R) i5-10210U CPU @ 1.60GHz
and 16GB of RAM, and a PC with 16 cores 11th Gen Intel(R) Core(TM) i7-11700 @
2.50GHz and 32 GB of RAM.

Baseline SVC algorithms have been measured on the notebook, while GWL times
have been determined with the PC.

##  5 Conclusion

Results confirm the theory, they’re not worse than the particular case. In
particular, it looks like that the preferred generalization scheme is the one
that gives more importance to patterns that are isolated, on 3 data-sets from
4 for the SVC case. Nevertheless it should be considered the fact concerning
the unique geometry of each data-set, so each generalization scheme should be
tested. The next step would be to leverage this method in order to learn the
weights. Perhaps this generalization could be employed in other contests such
as SVR, multi-class classification, and other MLPs.

## References

  * [1] ”Zhao Sen et al.”, _Metric-Optimized Example Weights_ , 2019.
  * [2] Aiolli F, Sperduti A., _An efficient SMO-like algorithm for multiclass SVM_ , Proceedings of the 12th IEEE Workshop on Neural Networks for Signal Processing, pp, 297–306, 2002/9/6
  * [3] Book: _Deep Learning_ , MIT Press, 2016
  * [4] Book: _Statistical Learning Theory_ WILEY, 1998
  * [5] Portera F., _A generalized quadratic loss for SVM and Deep Neural Networks_ LOD 2020 Conference work