Joint queue length distribution of multi-class, single-server queues with preemptive priorities

Abstract

In this paper we analyze an M/M/1 queueing system with an arbitrary number of customer classes, with class-dependent exponential service rates and preemptive priorities between classes. The queuing system can be described by a multi-dimensional Markov process, where the coordinates keep track of the number of customers of each class in the system. Based on matrix-analytic techniques and probabilistic arguments, we develop a recursive method for the exact determination of the equilibrium joint queue length distribution. The method is applied to a spare parts logistics problem to illustrate the effect of setting repair priorities on the performance of the system. We conclude by briefly indicating how the method can be extended to an M/M/1 queueing system with non-preemptive priorities between customer classes.

1 Introduction

We consider a single-server queueing system shared by N customer classes, numbered $1,\ldots ,N$. The class index n indicates the priority rank; class 1 has the lowest priority and class N has the highest priority. The arrival process of class-n customers is a Poisson process with rate $\lambda _n$. The service time of class-n customers is exponentially distributed with rate $\mu _n$. This system can be described by a multi-dimensional Markov process on the state space $\mathbb {N}_0^N$, where the coordinates keep track of the number of customers of each class in the system.

In this paper we present an exact method, based on matrix-analytic techniques [22, 23] to determine the equilibrium joint queue length distribution. In particular, it appears to be possible to avoid the use of infinite series and truncation of the state space. The crucial observation is that the Markov process, embedded on states in which there are no customers of priority classes higher than n, is of the M/G/1 type, where the number of class-n customers represents the class-n level. This is due to the fact that during excursions of the Markov process in which higher priority customers are present, any number of lower priority customers may arrive. Thus, a natural way to find the equilibrium joint queue length distribution is by recursive application of the theory of M/G/1-type Markov processes.

The joint queue length distribution is required in applications in the area of spare parts logistics and production. Specifically, our interest in the M/M/1 priority system with N classes arose from a spare parts logistics problem, where the joint queue length distribution is necessary for an exact performance analysis. This problem is discussed in Sect. 3.

Priority queueing systems have a long history (cf. [7, 8, 16]) and single- and multi-server priority queues received much attention. Most of the earlier studies concentrate on the transforms of marginal system characteristics such as the queue length and waiting time of a specific priority class. The focus on marginal system characteristics is also seen in recent work in [13, 28], where the domain of priority queueing systems with general arrival and service time distributions is treated.

Joint queue length distributions have first been studied in [19] using the matrix-geometric method [21] for an M/M/1 priority queueing system with two classes. This study spurred the observation made in [3, 29, 30] that the matrix-geometric method is a natural choice for studying priority queueing systems with a quasi-birth–death (QBD) structure. In these papers, the matrix-geometric method is generalized to systems with two priority classes, a Markovian arrival process and a phase-type service time distribution. In [4] the same matrix-geometric method is applied to a discrete-time N-class system, leading to an approximation of the joint equilibrium distribution, as the rate matrix R needs to be truncated for actual computation. An M/PH/1 non-preemptive priority system with N classes with different service rates per class is studied in [15], where an algorithm is derived using matrix-geometric techniques for the computation of the joint queue length distribution for three aggregated classes. The observation that is not made in [4, 15] is that lower priority customers see the queueing system as an M/G/1-type system, i.e., an M/M/1 system with an unreliable server (or vacations), where down times correspond to high-priority service interruptions. This observation is made and implemented in [12, 31], where the distribution of the down times are approximated by phase-type distributions, the first three moments of which are matched to the moments of high-priority service interruptions. However, only marginal queue length distributions are obtained.

There is also a number of papers studying the joint queue length distribution using alternative approaches. Generating functions are used in [9, 10] for the analysis of M/M/c priority queueing systems with two classes. Generating functions are also used in [20] for an M/M/c preemptive priority system with more than two classes. Here, customers of higher priority are aggregated, leading to an approximation of the equilibrium distribution. Later, [25–27] use a mixture of the matrix-geometric method and generating function technique to analyze preemptive and non-preemptive priority M/M/c queueing systems with two classes, where each class can have different types of customers. The mixture of the two methods leads to an approximation of the joint equilibrium distribution as the number of matrix operations has to be finite for actual computation.

The area of priority queueing systems still is an active field of research. More recently, priority queueing systems with impatient high-priority customers have been analyzed using generating functions [5]; by identifying simple Markov processes [6]; using a level-crossing method [14]; or using Laplace–Stieltjes transforms [17]. These systems have applications in, for example, telecommunication systems where voice messages need to be delivered quickly and have priority over data packets. An alternative to impatient customers are queueing systems where customers can reduce their sojourn time by transferring to a higher priority class. This allows impatient customers to be served earlier. In [32], bounds on the equilibrium distribution are given. The study of a queueing system with transferring customers is motivated by the potential application in the design of emergency departments. Here, patients are categorized in classes of different priority, where patients can transfer from a lower priority class to a higher priority class. Approximations for the first and second moment of the waiting time in an M/G/c non-preemptive priority queueing system with an arbitrary number of priority classes are given in [2].

Our main contribution is that we describe a method for the exact determination of the joint queue length distribution for a preemptive priority queueing system with an arbitrary number of classes and class-dependent service rates. We use the property that the embedded Markov process is of the M/G/1 type. A key to the approach is identifying first-passage probabilities which are computed by one-step analysis. We then recursively apply matrix-analytic methods related to M/G/1-type Markov processes and avoid the use of infinite series.

The remainder of the paper is organized as follows. In Sect. 2 we describe how the matrix-analytic method is applied to an N-class preemptive priority single-server system. To ease the understanding of the method in general and highlight the recursive nature, we first treat the two- and three-class systems in Sects. 2.1 and 2.2, respectively. Next, in Sect. 3 we present the application in spare parts logistics where the joint queue length distribution is needed for an exact analysis. In the final section we conclude by indicating how to extend the method to non-preemptive priority rules.

2 Matrix-analytic method

The M/M/1 preemptive priority system can be described by a Markov process with states $(q_N,\ldots ,q_1)$, where $q_n$ denotes the number of class-n customers in the system. State transitions are triggered by arrival and service completions. Class-n customers enter at an exponential rate $\lambda _n$, triggering a transition from $(q_N,\ldots ,q_1)$ to state $(q_N,\ldots ,q_n+1,\ldots ,q_1)$, and if $q_N = \cdots = q_{n+1} = 0$ and $q_n > 0$, class-n customers are served at an exponential rate $\mu _n$, which leads to a transition from $(0,\ldots ,0,q_n,\ldots ,q_1)$ to $(0,\ldots ,0,q_n-1,\ldots ,q_1)$. Throughout the paper we assume that the system is stable, i.e., the traffic intensity $\rho$ is less than 1 (see, for example, [11]):

$$\begin{aligned} \rho := \sum _{i = n}^N \lambda _n / \mu _n < 1, \end{aligned}$$

(1)

and we denote by $p(q_N,\ldots ,q_1)$ the equilibrium probability of being in state $(q_N,\ldots ,q_1)$. To ease notation, let us introduce $\lambda := \sum _{n = 1}^N \lambda _n$. We propose to use the matrix-analytic method for $M{\!/\!}G{\!/\!}1$ structured systems to exactly and recursively calculate the joint queue length probabilities $p(q_N,\ldots ,q_1)$, starting from $p(0,\ldots ,0) = 1 - \rho$. Key to this approach are first-passage probabilities, that can be determined through one-step analysis. In fact, the first-passage probabilities are the elements of the auxiliary matrix G of the matrix-analytic method. However, rather than determining the infinite matrix G using matrix equations, we recursively determine its elements using scalar equations, derived by exploiting the skip-free property of this Markov process. To highlight the recursive nature of the method, we first treat the two- and three-class systems.

Fig. 1

Transition rate diagrams of the two-class system. a Transition rate diagram of the two-class system, b Embedded on class-2 level $q_2 = 0$

2.1 Two-class system

The transition rate diagram of the two-class system depicted in Fig. 1a shows the two-class system is a QBD process with class-2 levels $q_2$ defined as the set of states with $q_2$ high-priority customers. To calculate the probabilities $p(q_2,q_1)$, we propose to exploit the $M{\!/\!}G{\!/\!}1$ structure of this Markov process, instead of its $G\!/\!M\!/\!1$ structure as done by [19]. Instrumental in the calculation of $p(q_2,q_1)$ are the first-passage probabilities $g_{2;i_1}$, instead of the elements of the rate matrix as in [19]. The first-passage probability $g_{2;i_1}$ is defined as the probability that, starting at class-2 level $q_2 > 0$ in state $(q_2,q_1)$, the first passage to class-2 level $q_2 - 1$ happens in state $(q_2-1,q_1 + i_1)$. Note that $g_{2;i_1}$ does not depend on the starting state $(q_2,q_1)$ and can be interpreted as the probability that $i_1$ class-1 customers arrive during a busy period of class-2 customers. By one-step analysis we get

$$\begin{aligned}&\mu _2 - (\lambda + \mu _2)g_{2;0} + \lambda _2 g_{2;0}^2 = 0, \quad i_1 = 0, \end{aligned}$$

(2)

$$\begin{aligned}&- (\lambda + \mu _2)g_{2;i_1} + \lambda _1 g_{2;i_1 - 1} + \lambda _2 \sum _{j_1 = 0}^{i_1} g_{2;j_1} g_{2;i_1-j_1} = 0, \quad i_1 > 0. \end{aligned}$$

(3)

So $g_{2;i_1}$ can be recursively calculated, starting from $g_{2;0}$, which follows from (2),

$$\begin{aligned} g_{2;0} = \frac{1}{2\lambda _2} \Bigl ( \lambda + \mu _2 - \bigl ((\lambda + \mu _2)^2 - 4\lambda _2\mu _2 \bigr )^{\frac{1}{2}} \Bigr ). \end{aligned}$$

(4)

To calculate $p(q_2,q_1)$, we use the following equation for excursions starting at class-2 level $q_2$ to levels higher than $q_2$ ending at first return to class-2 level $q_2$. The number of excursions per time unit that end in state $(q_2,q_1)$ is equal to $p(q_2+1,q_1)\mu _2$, but this number is also equal to the excursions starting from class-2 level $q_2$ per time unit that end in state $(q_2,q_1)$. The number of excursions per time unit that start in state $(q_2,q_1-i_1)$ is equal to $p(q_2,q_1-i_1)\lambda _2$, a fraction $g_{2;i_1}$ of which end in $(q_2,q_1)$. Hence,

$$\begin{aligned} p(q_2+1,q_1)\mu _2 = \sum _{i_1 = 0}^{q_1} p(q_2,q_1 - i_1) \lambda _2 g_{2;i_1}, \quad q_2,q_1 \ge 0, \end{aligned}$$

(5)

from which all probabilities can be recursively calculated, once the boundary probabilities $p(0,q_1)$ are known. The probabilities $p(0,q_1)$ can be determined by considering the Markov process embedded on class-2 level 0. The transition rate diagram of the embedded Markov process is shown in Fig. 1b. Note that the embedded Markov process has an $M{\!/\!}G{\!/\!}1$ structure with class-1 levels $q_1$ defined as the set of states with $q_1$ class-1 customers (and no class-2 customers). To formulate the analog of (5), we introduce $f_{2;i_1}$, which is the probability that, starting in state $(1,q_1)$, the first passage to class-1 levels less than or equal to $q_1 + i_1$ happens in state $(0,q_1 + i_1)$. In this case, this first-passage probability is equal to the probability that during a busy period of class-2 customers, at least $i_1$ class-1 customers arrive. So

$$\begin{aligned} f_{2;i_1} = 1 - \sum _{j_1 = 0}^{i_1 - 1} g_{2;j_1}. \end{aligned}$$

(6)

Then, similar to (5), we have

$$\begin{aligned} p(0,q_1 + 1)\mu _1 = p(0,q_1)\lambda _1 + \sum _{i_1 = 0}^{q_1} p(0,q_1 - i_1) \lambda _2 f_{2;i_1+1}, \quad q_1 \ge 0, \end{aligned}$$

(7)

which can be used to calculate all boundary probabilities, starting from the probability of an empty system $p(0,0) = 1 - \rho$.

Fig. 2

Transition rate diagram of the three-class system. a $q_3 > 0$, b $q_3 = 0$

2.2 Three-class system

The transition rate diagram of the three-class system is shown in Fig. 2a and b. This system can be described by a QBD process with class-3 levels $q_3$ defined as the set of states with $q_3$ high-priority customers. Let $g_{3;i_2,i_1}$ be the probability that, starting at class-3 level $q_3 > 0$ in state $(q_3,q_2,q_1)$, the first passage to class-3 level $q_3 - 1$ happens in state $(q_3 - 1,q_2 + i_2, q_1 + i_1)$. Note that $g_{3;i_2,i_1}$ can be interpreted as the probability that $i_2$ class-2 and $i_1$ class-1 customers arrive during a busy period of high-priority class-3 customers. By one-step analysis,

$$\begin{aligned}&\mu _3 - ( \lambda + \mu _3 ) g_{3;0,0} + \lambda _3 g_{3;0,0}^2 = 0, \quad i_2,i_1 = 0, \end{aligned}$$

(8)

$$\begin{aligned}&\quad - ( \lambda + \mu _3 )g_{3;i_2,i_1} + \lambda _1 g_{3;i_2,i_1-1} + \lambda _2 g_{3;i_2-1,i_1} \nonumber \\&\quad + \lambda _3 \sum _{j_2 = 0}^{i_2} \sum _{j_1 = 0}^{i_1} g_{3;j_2,j_1} g_{3;i_2-j_2,i_1-j_1} = 0, \quad i_2 + i_1 > 0, \end{aligned}$$

(9)

where, by convention, $g_{3;i_2,i_1} = 0$ if $i_2 < 0$ or $i_1 < 0$. From (9) the probabilities $g_{3;i_2,i_1}$ can be recursively calculated, starting from $g_{3;0,0}$, which follows from (8),

$$\begin{aligned} g_{3;0,0} = \frac{1}{2\lambda _3} \Bigl ( \lambda + \mu _3 + \bigl ( ( \lambda + \mu _3 )^2 - 4 \lambda _3 \mu _3 \bigr )^{\frac{1}{2}} \Bigr ). \end{aligned}$$

(10)

Similar to (5), we have

$$\begin{aligned}&p(q_3+1,q_2,q_1)\mu _3 \nonumber \\&\quad = \sum _{i_2 = 0}^{q_2} \sum _{i_1 = 0}^{q_1} p(q_3,q_2-j_2,q_1-j_1) \lambda _3 g_{3;i_2,i_1}, \quad q_3,q_2,q_1 \ge 0, \end{aligned}$$

(11)

which can be utilized to calculate all probabilities, once the boundary probabilities $p(0,q_2,q_1)$ are known.

Fig. 3

Transition rate diagram of two embedded Markov processes of the three-class system. a Embedded on class-3 level 0, b Embedded on the axis $q_3 = q_2 = 0$

To determine $p(0,q_2,q_1)$ we proceed by considering the Markov process embedded on class-3 level 0, which is of the $M{\!/\!}G{\!/\!}1$ type, with class-2 levels $q_2$ defined as the set of states with $q_2$ class-2 customers (and no class-3 customers). Its transition rate diagram is depicted in Fig. 3a. The first-passage probabilities $g_{2;i_1}$ for the embedded Markov process are defined as the probability that, when starting at class-2 level $q_2 > 0$ in state $(0,q_2,q_1)$, the first passage to class-2 level $q_2 - 1$ happens in state $(0,q_2-1,q_1+i_1)$. Further, the first-passage probabilities $g_{3;i_1}$ are defined as the probability that, when starting in state $(1,q_2-1,q_1)$, the first passage to class-2 level $q_2 - 1$ happens in state $(0,q_2-1,q_1+i_1)$. Observe that $g_{k;i_1}$ is the probability that $i_1$ class-1 customers arrive during a busy period of higher priority (class-2 and class-3) customers, that starts with the arrival of a class-k customer, for $k = 2,3$. Notice the difference between the first passage probabilities $g_{3;i_1}$ and $g_{3;i_2,i_1}$. The number of indices after the semicolon in the subscript is related to what level the Markov process is embedded on, as illustrated in Fig. 3. By one-step analysis we get for $k = 2,3$,

$$\begin{aligned}&\mu _k - ( \lambda + \mu _k ) g_{k;0} + \sum _{m = 2}^3 \lambda _m g_{m;0} g_{k;0} = 0, \quad i_1 = 0, \end{aligned}$$

(12)

$$\begin{aligned}&- ( \lambda + \mu _k ) g_{k;i_1} + \lambda _1 g_{k;i_1-1} + \sum _{m = 2}^3 \lambda _m \sum _{j_1 = 0}^{i_1} g_{m;j_1} g_{k;i_1-j_1} = 0, \quad i_1 > 0. \end{aligned}$$

(13)

From equations (13), both $g_{2;i_1}$ and $g_{3;i_1}$ can be recursively calculated, with $g_{2;0}$ and $g_{3;0}$ being the minimal non-negative solution of (12). To solve (12) we introduce $B_{k}$, which is the Laplace–Stieltjes transform (LST) of the service time of a class-k customer, and $BP_{k}$, the LST of a high-priority (class-2 and class-3) busy period initiated by a class-k customer. Then $g_{2;0}$ and $g_{3;0}$ can be calculated from, see [18, Sect. 5.8],

$$\begin{aligned} g_{k;0}&= BP_{k}(\lambda _1) = B_{k}(\lambda _1 + (\lambda _2 + \lambda _3)(1-BP_{2,3}(\lambda _1))) \nonumber \\&= \frac{\mu _k}{\mu _k + \lambda _1 + (\lambda _2 + \lambda _3)(1-BP_{2,3}(\lambda _1))}, \quad k = 2,3, \end{aligned}$$

(14)

where $BP_{2,3}(s)$ is the LST of a high-priority busy period initiated by a class-2 or a class-3 customer, which is equal to the LST of the busy period in an $M\!/\!H_2/\!1$ queue with class-2,3 customers,

$$\begin{aligned} BP_{2,3}(s) = \sum _{m = 2}^3 \frac{\lambda _m}{\lambda _2 + \lambda _3} \frac{\mu _m}{\mu _m + s + (\lambda _2 + \lambda _3)(1-BP_{2,3}(s))}, \quad s \ge 0. \end{aligned}$$

(15)

To formulate the analog of (11) for $p(0,q_2,q_1)$, we introduce the first-passage probabilities $f_{3;i_2,i_1}$ defined as the probability that, when starting in state $(1,q_2,q_1)$, the first passage to class-2 levels less than or equal to $q_2 + i_2$ happens in state $(0,q_2+i_2,q_1+i_1)$. The probability $f_{3;i_2,i_1}$ can be interpreted as the probability that at the end of a busy period of class-3 customers, there have been at least $i_2$ class-2 arrivals, and then, when the server brings down this number to $i_2$, the total number of class-1 arrivals (from the start of the busy period of class-3 customers) has been $i_1$. Hence, we can express $f_{3;i_2,i_1}$ as an infinite sum,

$$\begin{aligned} f_{3;i_2,i_1} = \sum _{m = 0}^\infty \sum _{\begin{array}{c} j_0,\ldots ,j_m \ge 0 \\ j_0 + \cdots + j_m = i_1 \end{array}} g_{3;i_2+m,j_0} g_{2;j_1} \cdots g_{2;j_m}. \end{aligned}$$

(16)

Before elaborating on the computation of $f_{3;i_2,i_1}$, we proceed to derive an equation for $p(0,q_2,q_1)$ by considering excursions to class-2 levels higher than $q_2$ that start at class-2 level $q_2$ or lower, and end at first return to class-2 level $q_2$ in state $(0,q_2,q_1)$. The number of excursions per time unit that end in state $(0,q_2,q_1)$ is equal to $p(0,q_2+1,q_1)\mu _2$. This number is also equal to the excursions starting from class-2 level $q_2$ or lower per time unit that end in state $(0,q_2,q_1)$. A fraction $g_{2;i_1}$ of the excursions starting in $(0,q_2,q_1-i_1)$ by a class-2 arrival end in $(0,q_2,q_1)$. Excursions to class-2 levels higher than $q_2$ starting in state $(0,q_2-i_2,q_1-i_1)$ by a class-3 arrival reach, with probability $f_{3;i_2,i_1} - g_{3;i_2,i_1}$, class-2 level $q_2$ in state $(0,q_2,q_1)$ at first return to class-2 level $q_2$. Note that $g_{3;i_2,i_1}$ needs to be subtracted, since with probability $g_{3;i_2,i_1}$ class-2 level $q_2$ is reached but not yet exceeded. Hence,

$$\begin{aligned}&p(0,q_2+1,q_1)\mu _2 = \sum _{i_1 = 0}^{q_1} p(0,q_2,q_1-i_1) \lambda _2 g_{2;i_1} \nonumber \\&\quad + \sum _{i_2 = 0}^{q_2} \sum _{i_1 = 0}^{q_1} p(0,q_2-i_2,q_1-i_1) \lambda _3 (f_{3;i_2,i_1} - g_{3;i_2,i_1}), \quad q_2,q_1 \ge 0, \end{aligned}$$

(17)

from which $p(0,q_2,q_1)$ can be recursively calculated, once the boundary probabilities $p(0,0,q_1)$ are known. To determine $p(0,0,q_1)$ we consider the Markov process embedded on the axis $q_3 = q_2 = 0$, the transition rate diagram of which is depicted in Fig. 3b, with class-1 levels $q_1$ defined as the set of states with $q_1$ class-1 customers (and no class-2 or class-3 customers). To finally formulate the equations for $p(0,0,q_1)$ we define $f_{k;i_1}$ as the probability that, when starting in state $(0,1,q_1)$ if $k = 2$ and starting in state $(1,0,q_1)$ if $k = 3$, the first passage to class-1 levels less than or equal to $q_1 + i_1$ happens in state $(0,0,q_1 + i_1)$. Similar to the two-class system, this first-passage probability is equal to the probability that at least $i_1$ class-1 customers arrive during a busy period of class-2, 3 customers, initiated by a class-k customer. So, for $k = 2,3$,

$$\begin{aligned} f_{k;i_1} = 1 - \sum _{j_1 = 0}^{i_1 - 1} g_{k;j_1}. \end{aligned}$$

(18)

Then, similar to (7), we have

$$\begin{aligned}&p(0,0,q_1+1)\mu _1 = p(0,0,q_1)\lambda _1 \nonumber \\&\quad + \sum _{i_1 = 0}^{q_1} p(0,0,q_1-i_1)(\lambda _2 f_{2;i_1+1} + \lambda _3 f_{3;i_1+1}), \quad q_1 \ge 0. \end{aligned}$$

(19)

This equation can be used to recursively calculate $p(0,0,q_1)$, with initially $p(0,0,0) = 1 - \rho$.

We now turn to the calculation of the first-passage probabilities $f_{3;i_2,i_1}$. To avoid evaluation of the infinite sums in (16), we again employ one-step analysis, yielding for $i_2 > 0$ and $i_1 \ge 0$,

$$\begin{aligned}&-( \lambda + \mu _3 ) f_{3;i_2,i_1} + \lambda _1 f_{3;i_2,i_1-1} + \lambda _2 f_{3;i_2-1,i_1} \nonumber \\&\quad + \lambda _3 \Bigl ( \sum _{j_2 = 0}^{i_2-1} \sum _{j_1 = 0}^{i_1} g_{3;j_2,j_1} f_{3;i_2-j_2,i_1-j_1} + \sum _{j_1 = 0}^{i_1} f_{3;i_2,j_1} g_{3;i_1-j_1} \Bigr ) = 0, \end{aligned}$$

(20)

where, by convention, $f_{3;i_2,i_1} = 0$ if $i_1 < 0$. The first-passage probabilities $f_{3;i_2,i_1}$ can be recursively calculated using the equations (20), starting with $f_{3;0,i_1} = g_{3;i_1}$. The last two terms in (20) need some explanation: this is the probability of first passage to class-2 levels less than or equal to $q_2 + i_2$ in state $(0,q_2+i_2,q_1+i_1)$ when starting an excursion in state $(2,q_2,q_1)$, so with two instead of one class-3 customer. Now imagine that the second class-3 customer enters service when the busy period generated by the first class-3 customer finishes. The first term corresponds to the event that the number of class-2 arrivals during the busy period generated by the first class-3 customer is $j_2 < i_2$, so that the number of class-2 arrivals during the second busy period should be at least $i_2 - j_2$. The second term corresponds to the event that the number of class-2 arrivals during the first busy period is $j_2 \ge i_2$. The surplus number $j_2 - i_2$ of class-2 customers should be served after the busy period generated by the second class-3 customer. The duration of the excursion will not be altered if these class-2 customers enter service (as well as any higher priority customer arriving during their service) before the second class-3 customer. Then $f_{3;i_2,j_1}$ is the probability that the number of class-1 arrivals is $j_1$ when the last surplus class-2 customer completes service, and thus, the number of class-1 arrivals during the busy period generated by the second class-3 customer should be exactly equal to $i_1 - j_1$. Note that the busy period generated by this second class-3 customer includes class-3 and class-2 customers, since each arriving class-2 customer is surplus. So the probability of exactly $i_1 - j_1$ class-1 arrivals is $g_{3;i_1-j_1}$.

2.3 N-class system

We now extend the approach for obtaining the stationary distribution of the three-class system to an N-class system. Since we are dealing with N classes, we need some accommodating notation. We introduce ${\mathbf {i}}^{(n)} = (i_n,i_{n-1},\ldots ,i_1)$, ${\mathbf {q}}^{(n)} = (q_n,q_{n-1},\ldots ,q_1)$, ${\mathbf {j}}^{(n)}$ is vector-index of length n, ${\mathbf {0}}^{(n)}$ is the zero vector of length n, and ${\mathbf {e}}^{(n)}_{k}$ denotes a vector of zeros of length n with a 1 at position $n + 1 - k$. Class-n level $q_n$ denotes the set of states with $q_n$ class-n customers and no customers of higher classes.

Once again, we have two types of first-passage probabilities. The first type is the first-passage probability $g_{k;{\mathbf {i}}^{(n)}}, ~ k \ge n + 1$ defined as the probability that, when starting in state $(0,\ldots ,0,q_{n+1}-1,{\mathbf {q}}^{(n)})+{\mathbf {e}}^{(N)}_{k}$, the first passage to class-$(n + 1)$ level $q_{n+1} - 1$ happens in state $(0,\ldots ,0,q_{n+1}-1,{\mathbf {q}}^{(n)} + {\mathbf {i}}^{(n)})$. Second, $f_{k;{\mathbf {i}}^{(n)}}, ~ k \ge n + 1$ is the probability that, when starting in state $(0,\ldots ,0,{\mathbf {q}}^{(n)}) + {\mathbf {e}}^{(N)}_{k}$, the first passage to class-n levels less than or equal to $q_{n}+i_{n}$ happens in state $(0,\ldots ,0,{\mathbf {q}}^{(n)}+{\mathbf {i}}^{(n)})$.

We first describe how to obtain the first-passage probabilities, followed by the computation of the equilibrium probabilities. Note that $f_{k;0,{\mathbf {i}}^{(n-1)}} = g_{k;{\mathbf {i}}^{(n-1)}}$. By one-step analysis we get, for $n = N-1,N-2,\ldots ,1$ and $k \ge n + 1$,

$$\begin{aligned}&\mu _k - (\lambda + \mu _k)g_{k;{\mathbf {0}}^{(n)}} + \sum _{m = n + 1}^N \lambda _m g_{m;{\mathbf {0}}^{(n)}} g_{k;{\mathbf {0}}^{(n)}} = 0, \end{aligned}$$

(21)

$$\begin{aligned}&\quad -(\lambda + \mu _k)g_{k;{\mathbf {i}}^{(n)}} + \sum _{m = 1}^{n} \lambda _m g_{k;{\mathbf {i}}^{(n)}-{\mathbf {e}}^{(n)}_{m}} \nonumber \\&\quad + \sum _{m = n + 1}^N \lambda _m \sum _{{\mathbf {j}}^{(n)} = {\mathbf {0}}^{(n)}}^{{\mathbf {i}}^{(n)}} g_{m;{\mathbf {j}}^{(n)}} g_{k;{\mathbf {i}}^{(n)}-{\mathbf {j}}^{(n)}} = 0, \quad {\mathbf {i}}^{(n)} > {\mathbf {0}}^{(n)}. \end{aligned}$$

(22)

From (22), all $g_{k;{\mathbf {i}}^{(n)}}$ with $k \ge n + 1$ and n fixed can be calculated, with $g_{k;{\mathbf {0}}^{(n)}}$ computed as

$$\begin{aligned} g_{k;{\mathbf {0}}^{(n)}}&= BP_{k}\Bigl (\sum _{m = 1}^{n} \lambda _m \Bigr ) \nonumber \\&= B_{k}\Bigl (\sum _{m = 1}^{n} \lambda _m + \sum _{m = n+1}^N \lambda _m (1 - BP_{n+1,\ldots ,N}(\sum _{m = 1}^{n} \lambda _m )) \Bigr ) \nonumber \\&= \frac{\mu _k}{\mu _k + \sum _{m = 1}^{n} \lambda _m + \sum _{m = n+1}^N \lambda _m (1 - BP_{n+1,\ldots ,N}(\sum _{m = 1}^{n} \lambda _m))}, \end{aligned}$$

(23)

where $BP_{n+1,\ldots ,N}(s)$ is the LST of a high-priority (class-$n+1$ and higher) busy period, which is equal to the LST of the busy period in an $M\!/\!H_{N-n}/1$ queue with class-$n+1,\ldots ,N$ customers,

$$\begin{aligned} BP_{n+1,\ldots ,N}(s)&= \sum _{m = n+1}^{N} \frac{\lambda _m}{\sum _{l = n+1}^N \lambda _l}\nonumber \\&\quad \cdot \frac{\mu _m}{\mu _m + s + \sum _{l = n+1}^N \lambda _l(1-BP_{n+1,\ldots ,N}(s))}, \quad s \ge 0. \end{aligned}$$

(24)

The first-passage probabilities $f_{k;{\mathbf {i}}^{(n)}}$ with $n = N-1,\ldots ,2$ and $k \ge n + 1$ follow from one-step analysis similar to (20), with $i_n > 0$,

$$\begin{aligned}&-(\lambda + \mu _k)f_{k;{\mathbf {i}}^{(n)}} + \sum _{m = 1}^{n} \lambda _m f_{k;{\mathbf {i}}^{(n)}-{\mathbf {e}}^{(n)}_{m}} \nonumber \\&\quad +\sum _{m = n+1}^N \lambda _m \Bigl ( \sum _{j_{n} = 0}^{i_{n} - 1} \sum _{{\mathbf {j}}^{(n-1)} = {\mathbf {0}}^{(n-1)}}^{{\mathbf {i}}^{(n-1)}} g_{m;j_n,{\mathbf {j}}^{(n-1)}} f_{k;i_n-j_n} \nonumber \\&\quad + \sum _{{\mathbf {j}}^{(n-1)} = {\mathbf {0}}^{(n-1)}}^{{\mathbf {i}}^{(n-1)}} f_{m;i_{n},{\mathbf {j}}^{(n-1)}} g_{k;{\mathbf {i}}^{(n-1)}-{\mathbf {j}}^{(n-1)}} \Bigr ) = 0, \quad {\mathbf {i}}^{(n-1)} \ge {\mathbf {0}}^{(n-1)}. \end{aligned}$$

(25)

The last two terms in (25) describe the probability of first passage to class-n levels less than or equal to $q_{n} + i_{n}$ in state $({\mathbf {0}}^{(N-n)},{\mathbf {q}}^{(n)}+{\mathbf {i}}^{(n)})$ when starting an excursion in state $({\mathbf {0}}^{(N-n)},{\mathbf {q}}^{(n)}) + {\mathbf {e}}^{(N)}_{k} + {\mathbf {e}}^{(N)}_{m}$, so with one class-k and one class-m customer. Note that we act as if the class-k customer enters service when the high-priority busy period generated by the class-m customer finishes. This is feasible, since the order in which the customers are served does not alter the duration of a high-priority busy period, cf. (20). The remaining first-passage probabilities for the case $n = 1$ are computed as, for $k \ge 2$,

$$\begin{aligned} f_{k;i_1} = 1 - \sum _{j_1 = 0}^{i_1 - 1} g_{k;j_1}, \quad i_1 > 0. \end{aligned}$$

(26)

The equilibrium probabilities of the N-class system follow again by counting excursions as done for the two- and three-class systems. The number of excursions per time unit that end in state $({\mathbf {0}}^{(N-n)},{\mathbf {q}}^{(n)})$ is equal to $p({\mathbf {0}}^{(N-n)},q_n+1,{\mathbf {q}}^{(n-1)})\mu _n$. This number is also equal to the excursions starting from class-n level $q_n$ or lower per time unit that end in state $({\mathbf {0}}^{(N-n)},{\mathbf {q}}^{(n)})$. A fraction $g_{n;{\mathbf {i}}^{(n-1)}}$ of the excursions starting in $({\mathbf {0}}^{(N-n)},q_n,{\mathbf {q}}^{(n-1)}-{\mathbf {i}}^{(n-1)})$ by a class-n arrival end in $({\mathbf {0}}^{(N-n)},{\mathbf {q}}^{(n)})$. Excursions to class-n levels higher than $q_n$ starting in state $({\mathbf {0}}^{(N-n)},{\mathbf {q}}^{(n)}-{\mathbf {i}}^{(n)})$ by a class-$m, ~ m = n + 1,\ldots ,N$ arrival reach, with probability $f_{m;{\mathbf {i}}^{(n)}} - g_{m;{\mathbf {i}}^{(n)}}$, level $q_n$ in state $({\mathbf {0}}^{(N-n)},{\mathbf {q}}^{(n)})$ at first return to class-n level $q_n$. We have for $n = 1,2,\ldots ,N$,

$$\begin{aligned}&p({\mathbf {0}}^{(N-n)},q_n+1,{\mathbf {q}}^{(n-1)})\mu _n \nonumber \\&\quad = \sum _{{\mathbf {i}}^{(n)} = {\mathbf {0}}^{(n)}}^{{\mathbf {q}}^{(n)}} p({\mathbf {0}}^{(N-n)},{\mathbf {q}}^{(n)}-{\mathbf {i}}^{(n)}) \sum _{m = n+1}^N \lambda _m (f_{m;{\mathbf {i}}^{(n)}} - g_{m;{\mathbf {i}}^{(n)}}) \nonumber \\&\quad \quad + \sum _{{\mathbf {i}}^{(n-1)} = {\mathbf {0}}^{(n-1)}}^{{\mathbf {q}}^{(n-1)}} p({\mathbf {0}}^{(N-n)},q_n,{\mathbf {q}}^{(n-1)}-{\mathbf {i}}^{(n-1)}) \lambda _n g_{n;{\mathbf {i}}^{(n-1)}}, \quad {\mathbf {q}}^{(n)} \ge {\mathbf {0}}^{(n)}, \end{aligned}$$

(27)

which can be solved recursively, starting from $p({\mathbf {0}}^{(N)}) = 1 - \rho$. Note that for $n = 1$, the second term on the right-hand side of (27) becomes $p({\mathbf {0}}^{(N-1)},q_1)\lambda _1$ and for $n = N$, the first term on the right-hand side reduces to 0.

Remark 1

The above algorithm to determine the equilibrium probabilities involves subtractions in some equations, see, for example, (26), which may possibly lead to loss of significant digits and instability. However, in all experiments we observed numerically stable results.

3 Application in spare parts logistics

Our interest in the joint queue length distribution arose from a spare parts supply problem for repairable parts sharing the same repair shop. For this problem, we apply our method, based on the matrix-analytic approach, to demonstrate the influence of assigning repair priorities on the performance of the system.

There are M identical machines and each machine contains three different subsystems, numbered 1, 2, 3. Each subsystem n consists of $Z_n$ identical parts in parallel. We refer to the parts of subsystem n as parts of Stock-Keeping Unit n (SKU n). For each subsystem, $k_n < Z_n$ parts have to function, i.e., we have redundancy, and the redundant parts are in “cold standby.” This is called a “$k_n$-out-of-$Z_n$” setup. We have $k_n$ functioning parts per subsystem and only these parts are subject to failure. When one part fails, another one can immediately take over the necessary functions. An example of such a subsystem is the board computer of an airplane, where this critical component is duplicated and in an idle mode to accommodate possible failures; here, $k_n = 1$ and $Z_n = 2$. Other typical systems with this structure can be found in [24].

A machine is only working when all three subsystems are working. When one of the functioning parts fails, a redundant part takes over its function and a service engineer takes a new part from a stock of parts and replaces the failed one. The failed part is then sent to a single-server repair facility. Part and repair requests are served on a first come first serve basis. The repair time for a part of SKU n is exponentially distributed with rate $\mu _n$; the delivery and replacement times are small and can be neglected. We assume that failures of parts of SKU n occur according to a Poisson process with rate $\lambda _n$. This approximation, which is the only one needed, is valid when M, the total number of machines in the system, is large and when the fraction of working machines is high. After repair the broken parts are assumed to be as good as new, and they are put back to stock. The stock of SKU n at time instant $t = 0$ is denoted by $S_n$. We call the amount $S_n$ the basestock level for SKU n parts. The system is shown in Fig. 4.

Fig. 4

Example of a simple spare parts supply system

Let us define the system availability as the average fraction of working machines:

$$\begin{aligned} A(S_1,S_2,S_3) = \frac{1}{M} \sum _{m = 1}^M {\mathbb {P}}(\text {Machine }m\hbox { is working}). \end{aligned}$$

(28)

The number of backorders of SKU n parts is given by $(q_n - S_n)^+$, where $(x)^+ = \max (0,x)$ and $q_n$ is the number of SKU n parts in repair. Define $E_n$ as the number of “empty” spots in a given subsystem n of any of the M machines. Then, by conditioning on the number of parts of SKU n in repair we obtain, with $s \le Z_n$,

$$\begin{aligned} {\mathbb {P}}(E_n = s \mid q_n \text { rep.}) = \left\{ \begin{array}{cl} 1, &{} q_n - S_n < s, ~ s = 0, \\ 0, &{} q_n - S_n < s, ~ s > 0, \\ \left( {\begin{array}{c}Z_n\\ s\end{array}}\right) \left( {\begin{array}{c}Z_n(M-1)\\ q_n-S_n-s\end{array}}\right) / \left( {\begin{array}{c}M Z_n\\ q_n - S_n\end{array}}\right) ,&q_n - S_n \ge s. \end{array} \right. \end{aligned}$$

(29)

In terms of the joint queue length distribution, the system availability can be written as

$$\begin{aligned} A(S_1,S_2,S_3) = \sum _{q_3,q_2,q_1 \ge 0} \Bigl ( \prod _{n = 1}^3 \sum _{s = 0}^{Z_n-k_n} {\mathbb {P}}(E_n = s \mid q_n \text { rep.}) \Bigr ) p(q_3,q_2,q_1). \end{aligned}$$

(30)

The expression (30) determines the system availability much better than other approximations proposed in the literature, for example, the system availability defined in [24] only uses information on the mean number of backorders. The matrix-analytic method makes it possible to use the detailed distribution of the number of parts in repair. To demonstrate the approach we execute a set of experiments with the following parameters: $M = 100$, and $Z_n = 4$, $k_n = 2$, $\lambda _n = n/300$, and $\mu _n = (4-n)/\beta$ for $n = 1,2,3$, where $\beta$ is chosen such that $\sum _{n = 1}^3 \lambda _n / \mu _n = \rho$.

We wish to compute the joint queue length distribution such that the sum $\sum _{q_3,q_2,q_1} p(q_3,q_2,q_1) > 1 - \epsilon$ with $\epsilon$ a small positive number. We do this by computing the equilibrium probabilities of the states in a discrete three-dimensional cuboid $\mathcal {C}$ with states $\{0,\ldots ,c_3\} \times \{0,\ldots ,c_2\} \times \{0,\ldots ,c_1\}$. For the sake of clarity, we briefly introduce the marginal queue length distribution of class-n customers as $p_n(\cdot )$. We specify the construction of $\mathcal {C}$ in more detail. The bound $c_3$ is computed from the M/M/1 system with only class-3 customers, such that $\sum _{q_3 = 0}^{c_3} p_3(q_3) > 1 - \epsilon$, which leads to $c_3 = \lceil \frac{\log {\epsilon }}{\log {\lambda _3/\mu _3}} - 1 \rceil$. The bound $c_2$ is obtained through a priority system with class-3, 2 customers such that the sum of the marginal probabilities for class-2 customers is very close to 1. That is, $\sum _{q_{2} = 0}^{c_{2}} p_2(q_2) > 1 - \epsilon$. Conveniently, the marginal queue length distribution $p_2(\cdot )$ can be derived directly from the joint equilibrium probabilities $p(0,\cdot )$ of the priority queueing system with class-3, 2 customers via the relation $\lambda _2 p_2(q_2-1) = \mu _2 p(0,q_2)$. Thus, this allows us to estimate the bound $c_2$ without having to compute all joint equilibrium probabilities $p(q_3,q_2)$. The final bound $c_1$ can be found iteratively until $\sum _{q_3,q_2,q_1} p(q_3,q_2,q_1) > 1 - \epsilon$ or using the same method as for the bound $c_2$. Naturally, this method of constructing $\mathcal {C}$ extends to an arbitrary number of classes.

In Table 1 we list the system availability according to (30) for different utilization rates of the repair shop and different priority assignments. The basestock levels $S_n$ depend on the mean queue lengths, i.e., we set $S_n = \lfloor {\mathbb {E}}[Q_n] \rfloor , ~ n = 1,2,3$, where $Q_n$ is the queue length of SKU n parts. The algorithm for the 3-class system was executed using Java 8.0 on a PC with an Intel Core i7-3770 CPU and 16 GB RAM. The computation times mentioned in Table 1 depend on the number of states with significant probability mass, i.e., on the load of the system, the priority assignment and, naturally, the parameter value of $\epsilon$. For these experiments, we have selected $\epsilon = 10^{-6}$.

Table 1 System availability for different combinations of repair shop utilizations and priority assignments

Fig. 5

Transition rate diagram of the non-preemptive M/M/1 priority system of the states $(0,q_1,1)$ with $q_1 > 0$ and including state (0, 0, 0). The dashed arrows indicate a transition to a state with $s = 2$, i.e., a state with a class-2 customer in service

Table 1 shows that we have a fast numerical method to compute the availability for different priority assignments and particular choices of the basestock levels. This method can easily be exploited in a procedure to optimize the priority assignment and basestock level, for example, in order to maximize system availability under a given budget for spare parts (cf. [1] which considers a slightly different setting with equal repair rates for all SKUs).

4 Conclusion and extensions

We have developed for the M/M/1 preemptive priority system with N customer classes and class-dependent service rates a method for the exact determination of the joint equilibrium queue length distribution. This method is based on the matrix-analytic method as the embedded Markov processes are of the $M{\!/\!}G{\!/\!}1$ type. Key to this approach are first-passage probabilities, computed by one-step analysis.

We applied the exact solution method to a spare parts logistics problem where repairable parts share the same repair shop, and showed that this method produces accurate results in the order of seconds.

We next sketch how the method can be extended to an M/M/1 non-preemptive priority system. In the non-preemptive case one identifies the customer currently in service by adding another variable to the state description. For the two-class system, the state description becomes $(q_2,q_1,s)$, where $s \in \{1,2\}$ indicates the class of the customer in service and $s = 0$ indicates no customer in service. By defining class-2 level $q_2$ as the set of states with $q_2$ class-2 customers, one can again count the number of excursions per time unit that start from class-2 level $q_2$ and reach levels higher than $q_2$ to finally end at state $(q_2,q_1,2)$. The states with a class-1 customer in service can only be reached from the states $(q_2,q_1,1)$ or (0, 0, 0), and thus, the equilibrium probabilities of these states can be recursively determined for $q_2 > 0$ immediately from the boundary probabilities of class-2 level 0, see Fig. 5. One finds the equilibrium probabilities of class-2 level 0, starting from $p(0,0,0) = 1 - \rho$, by embedding the Markov process on class-2 level 0 and again counting excursions. Notice that the approach is very similar to the one for the preemptive case and only requires the computation of equilibrium probabilities of the states $(q_2,q_1,1)$ as an additional step.