Time-dependent analysis of an $M/M/c$ preemptive priority system with two priority classes

We analyze the time-dependent behavior of an $M/M/c$ priority queue having two customer classes, class-dependent service rates, and preemptive priority between classes. More particularly, we develop a method that determines the Laplace transforms of the transition functions when the system is initially empty. The Laplace transforms corresponding to states with at least $c$ high-priority customers are expressed explicitly in terms of the Laplace transforms corresponding to states with at most $c - 1$ high-priority customers. We then show how to compute the remaining Laplace transforms recursively, by making use of a variant of Ramaswami's formula from the theory of $M/G/1$-type Markov processes. While the primary focus of our work is on deriving Laplace transforms of transition functions, analogous results can be derived for the stationary distribution: these results seem to yield the most explicit expressions known to date.


Introduction
Priority models with multiple servers constitute an important class of queueing systems, having applications in areas as diverse as manufacturing, wireless communication and the service industry. Studies of these models date back to at least the 1950s (see, for example, Cobham [7], Davis [8], and Jaiswal [15,16]), yet many properties of these systems still do not appear to be well understood. Recent work addressing priority models includes Sleptchenko et al. [27] and Wang et al. [29]. We refer the reader to [29] for more specific examples of applications of priority queueing models.
Our contribution to this stream of the literature is an analysis of the timedependent behavior of a Markovian multi-server queue with two customer classes, class-dependent service rates and preemptive priority between classes. To the best of our knowledge, the joint stationary distribution of the M/M/1 2-class preemptive priority system was first studied in Miller [24], who makes use of matrix-geometric methods to study the joint stationary distribution of the number of high-and lowpriority customers in the system. More particularly, in [24] this queueing system is modeled as a quasi-birth-and-death (QBD) process having infinitely many levels, with each level containing infinitely many phases. Miller then shows how to recursively compute the elements of the rate matrix of this QBD process: Once enough elements of this rate matrix have been found, the joint stationary distribution can be approximated by appropriately truncating this matrix.
This single-server model is featured in many works that have recently appeared in the literature. In Sleptchenko et al. [27] an exact, recursive procedure is given for computing the joint stationary distribution of an M/M/1 preemptive priority queue that serves an arbitrary finite number of customer classes. The M/M/1 2-class priority model is briefly discussed in Katehakis et al. [20], where they explain how the successive lumping technique can be used to study M/M/1 2-class priority models when both customer classes experience the same service rate. Interesting asymptotic properties of the stationary distribution of the M/M/1 2-class preemptive priority model can be found in the work of Li and Zhao [23].
Multi-server preemptive priority systems with two customer classes have also received some attention in the literature. One of the earlier references allowing for different service requirements between customer classes is Gail et al. [13]; see also the references therein. In [13], the authors derive the generating function of the joint stationary probabilities by expressing it in terms of the stationary probabilities associated with states where there is no queue. A combination of a generating function approach and the matrix-geometric approach is used in Sleptchenko et al. [26] to compute the joint stationary distribution of an M/M/c 2-class preemptive priority queue. The M/P H/c queue with an arbitrary number of preemptive priority classes is studied in Harchol-Balter et al. [14] using a recursive dimensionality reduction technique that leads to an accurate approximation of the mean sojourn time per customer class. Furthermore, in Wang et al. [29] the authors present a procedure for finding, for an M/M/c 2-class priority model, the generating function of the distribution of the number of low-priority customers present in the system in stationarity.
Our work deviates from all of the above approaches, in that we construct a procedure for computing the Laplace transforms of the transition functions of the M/M/c 2-class preemptive priority model. Our method first makes use of a slight tweak of the clearing analysis on phases (CAP) method featured in Doroudi et al. [10], in that we show how CAP can be modified to study Laplace transforms of transition functions. The specific dynamics of our priority model allow us to take the analysis a few steps further, by showing each Laplace transform can be expressed explicitly in terms of transforms corresponding to states contained within a strip of states that is infinite in only one direction. Finally, we show how to compute these remaining transforms recursively, by making use of a slight modification of Ramaswami's formula [25]. While the focus of our work is on Laplace transforms of transition functions, analogous results can be derived for the stationary distribution of the M/M/c 2-class preemptive priority model as well. We are not aware of any studies that obtain explicit expressions for the Laplace transforms of the transition functions, or even the stationary distribution, as we do here; these results seem to yield the most explicit expressions known to date.
The Laplace transforms we derive can easily be numerically inverted to retrieve the transition functions with the help of the algorithms of Abate and Whitt [3,4] or den Iseger [9]. These transition functions can be used to study-as a function of time-key performance measures such as the mean number of customers of each priority class in the system; the mean total number of customers in the system; or the probability that an arriving customer has to wait in the queue. These time-dependent performance measures can, for example, be used to analyze and dimension priority systems when one is interested in the behavior of such systems over a finite time horizon. Using the equilibrium distribution as an approximation of the time-dependent behavior to dimension the system can result in either over-or underdimensioning since it neglects transient effects such as the initial number of customers in the system. Prior to our work, the time-dependent behavior of multi-server queues could not be analyzed using any methods from previous work; until now one would have to resort to simulation in order to study the system's time-dependent behavior. Having explicit expressions for the Laplace transforms of the transition functions greatly simplifies the computation of some performance measures. For instance, these transforms yield explicit expressions for the Laplace transforms of the distribution of the number of low-priority customers in the system at time t.
We now present some numerical examples of the time-dependent performance measures, where we will make use of the notation introduced in Sect. 2. In Fig. 1, we plot the mean number of low-priority customers in the system as a function of time and we also depict the equilibrium values. Similarly, in Fig. 2 we plot the time-dependent and equilibrium delay probabilities for each priority class. The Laplace transforms used to obtain Figs. 1 and 2 can be computed numerically using the approach discussed in Sect. 5.5: here we used an error tolerance of = 10 −8 . Once these transforms have been found, numerical inversion can be done via the Euler summation algorithm of [3] where we again used an error tolerance of 10 −8 . From Fig. 1 we can also informally derive the mixing times of each scenario. It seems that the mixing time vastly increases with an increase in the load. As expected, in Fig. 2 we see that the delay probability of a high-priority customer is much lower than the delay probability of a low-priority customer. Furthermore, as time passes, the delay probability of the high-priority customer tends to the delay probability in an M/M/c queue with only high-priority customers, which was verified to be correct. Finally, in Table 1 we show   We use the numerical implementation outlined in Sect. 5.5 with accuracies = 10 −8 . Parameter settings are ρ 1 = 1/3, ρ 2 = 1/2 and c varies the computation times of the algorithm, which was implemented in Matlab and run on a 64-bit desktop with an Intel Core i7-3770 processor. The computation time scales reasonably well with the number of servers, and therefore the algorithm can be used to evaluate any practical instance.
This paper is organized as follows. Section 2 describes both the M/M/c 2-class preemptive priority queueing system, as well as the two-dimensional Markov process used to model the dynamics of this system. In the same section, we introduce relevant notation and terminology and detail the outline of the approach. In Sects. 3-5, we describe this approach for calculating the Laplace transforms of the transition functions. We discuss the simplifications in the single-server case in Sect. 6. In Sect. 7, we summarize our contributions and comment on the derivation of the stationary distribution. The appendices provide supporting results on combinatorial identities and single-server queues used in deriving the expressions for the Laplace transforms.

Model description and outline of approach
We consider a queueing system consisting of c servers, where each server processes work at unit rate. This system serves customers from two different customer classes, referred to here as class-1 and class-2 customers. The class index indicates the priority rank, meaning that among the servers, class-2 customers have preemptive priority over class-1 customers in service. Recall that the term 'preemptive priority' means that whenever a class-2 customer arrives at the system, one of the servers currently serving a class-1 customer immediately drops that customer and begins serving the new class-2 arrival, and the dropped class-1 customer waits in the system until a server is again available to receive further processing, i.e., the priority rule is preemptive resume. Therefore, if there are currently i class-1 customers and j class-2 customers in the system, the number of class-2 customers in service is min(c, j), while the number of class-1 customers in service is max(min(i, c − j), 0).
Class-n customers arrive in a Poisson manner with rate λ n , n = 1, 2, and the Poisson arrival processes of the two populations are assumed to be independent. Each class-n arrival brings an exponentially distributed amount of work with rate μ n , independently of everything else. We denote the total arrival rate by λ := λ 1 +λ 2 , the load induced by class-n customers as ρ n := λ n /(cμ n ), and the load induced by both customer classes as ρ := ρ 1 + ρ 2 .
The dynamics of this queueing system can be described with a continuous-time Markov chain (CTMC). For each t ≥ 0, let X n (t) represent the number of class-n customers in the system at time t, and define X (t) := (X 1 (t), X 2 (t)). Then, X := {X (t)} t≥0 is a CTMC on the state space S = N 2 0 . Given any two distinct elements x, y ∈ S, the element q(x, y) of the transition rate matrix Q associated with X denotes the transition rate from state x to state y. The row sums of Q are 0, meaning for each x ∈ S, q(x, x) = − y =x q(x, y) =: −q(x), where q(x) represents the rate of each exponential sojourn time in state x. For our queueing system, the nonzero transition rates of Q are given by We further associate with the Markov process X the collection of transition functions { p x,y (·)} x,y∈S , where for each x, y ∈ S (with possibly x = y) the function p x,y : [0, ∞) → [0, 1] is defined as Each transition function p x,y (·) has a Laplace transform π x,y (·) that is well-defined on the subset of complex numbers C + := {α ∈ C : Re(α) > 0} as We restrict our interest to transition functions of X when X (0) = (0, 0) with probability one (w.p.1), and so we drop the first subscript on both transition functions and Laplace transforms, i.e., p x (t) := p (0,0),x (t) for each t ≥ 0 and π x (α) := π (0,0),x (α) for each α ∈ C + . Our goal is to derive efficient numerical methods for calculating each Laplace transform π x (α), x ∈ S. We often refer to the Laplace transform π x (α) associated with the state x as the Laplace transform for state x. We sometimes refer to upper level U 0 as the vertical boundary. The union of upper levels is called the interior of the state space. Finally, the union is called the horizontal boundary or the horizontal strip of boundary states. Figure 4 depicts these sets. The indicator function 1{A} equals 1 if A is true and 0 otherwise. Given an arbitrary CTMC Z , we let E z [ f (Z )] represent the expectation of a functional of Z , conditional on Z (0) = z, and P z (·) denotes the conditional probability associated with E z [·]. In our analysis, it should be clear from the context what is being conditioned on when we write P z (·) or E z [·].
We will also need to make use of hitting-time random variables. We define for each set A ⊂ S, as the first time X makes a transition into the set A (so note X (0) ∈ A does not imply τ A = 0) and τ x should be understood to mean τ {x} .

Notation for M/M/1 queues
Most of the formulas we derive contain quantities associated with an ordinary M/M/1 queue. Given an M/M/1 queueing system with arrival rate λ and service rate μ, let Q λ,μ (t) denote the total number of customers in the system at time t. Under the measure P n (·), which, in this case, represents conditioning on Q λ,μ (0) = n, let B λ,μ denote the busy period duration induced by these customers. Under P 1 (·), the Laplace-Stieltjes transform of B λ,μ is given by Recall that under P n (·), B λ,μ is equal in distribution to the sum of n i.i.d. copies of B λ,μ under the measure P 1 (·); see, for example, [28, p. 32]. Thus, for each integer We will also need to make use of the following quantities in Sects. 5 and 6. Suppose {Λ θ (t)} t≥0 is a homogeneous Poisson process with rate θ that is independent of (2.8) Lemma 7 of Appendix C shows that the w (λ,μ,θ) i (α) terms satisfy a recursion, and, in Lemma 8 of Appendix C, we give explicit expressions for these terms by solving the recursion.
The following quantities associated with M/M/1 queues will appear at many places of the analysis. To increase readability, we adopt the notation used in [10] and define the quantities and . (2.10) Further results for M/M/1 queues are presented in Appendix C.

Outline of our approach
Our approach for computing the Laplace transforms of the transition functions of X when X (0) = (0, 0) w.p.1 is divided into three parts: 1. For each integer i ≥ 0, we use a slight modification of the CAP method [10] to write each Laplace transform for each state in U i , i.e., π (i,c−1+ j) (α), j ≥ 1, in terms of the Laplace transforms π (k,c−1) (α), 0 ≤ k ≤ i, as well as additional coefficients {v k,l } i≥k≥l≥0 that satisfy a recursion. 2. In Sect. 4, we obtain an explicit expression for the coefficients {v k,l } k≥l≥0 . This in turn shows that for each i ≥ 0, each Laplace transform for each state in U i , i.e., In Sect. 5, we derive a recursion with which we can determine the Laplace transforms for the states in the horizontal boundary. Specifically, we derive a modification of Ramaswami's formula [25] to recursively compute the remaining Laplace transforms The techniques we use to derive this recursion are exactly the same as the techniques recently used in [17] to study block-structured Markov processes. Only the Ramaswami-like recursion is needed to compute all Laplace transforms: Once the values for the Laplace transforms of the states in the horizontal boundary are known, all other transforms can be stated explicitly without using additional recursions.

A slight modification of the CAP method
The following theorem is used in multiple ways throughout our analysis. It appears in [ or, equivalently,
Using Theorem 1 with A = U c 0 we obtain, for j ≥ 1, From the transition rate diagram in Fig. 3, we find that the expectation in (3.3) can be interpreted as an expectation associated with an M/M/1 queue having arrival rate λ 2 and service rate cμ 2 . Indeed, τ U c 0 is equal in distribution to the minimum of the busy period-initialized by one customer-of this M/M/1 queue and an exponential random variable with rate λ 1 that is independent of the queue. Alternatively, τ U c 0 can be thought of as being equal in distribution to the busy period duration of an M/M/1 clearing model, with arrival rate λ 2 , service rate cμ 2 , and clearings that occur in a Poisson manner with rate λ 1 . Applying Lemma 6 of Appendix C shows that (3.4)

Laplace transforms for states within the interior
We next develop a recursion for the Laplace transforms for the states within the interior. First, we express the transforms in upper level U i in terms of the transforms in upper level U i−1 and in state (i, c −1). Second, we use this result to express the transforms in upper level U i in terms of the transforms for the states (0, and some additional coefficients. Employing again Theorem 1, now with A = U c i , yields, for i, j ≥ 1, The expectation has the same interpretation as the expectation in (3.3), except now the M/M/1 queue starts with k customers at time 0. Using Lemma 6 of Appendix C, we obtain where, for j, k ≥ 1, Substituting (3.8) into (3.6) and simplifying yields a recursion. Specifically, for i ≥ 0 and j ≥ 1, where the quantities {v i, j } i≥ j≥0 satisfy the following recursive scheme: Throughout we follow the convention that all empty sums, such as 0 k=1 (·), represent the number zero.

Deriving an explicit expression for v i, j
Theorem 2 suggests that the Laplace transforms for the states within the interior can be computed recursively: c. Once steps 2a. and 2b. are completed, all Laplace transform values π (i+1,c−1+ j) (α), j ≥ 1, are known. Our next result, i.e., Theorem 3, shows that for each i ≥ 0, the {v i,l } i≥l≥0 terms can be expressed explicitly in terms of π (k,c−1) (α), 0 ≤ k ≤ i. If our goal is to only compute π (i,c−1+ j) (α) for some large i, then Theorem 3 allows us to avoid computing all intermediate {v k,l } i−1≥k≥l≥0 terms, which means we can avoid computing an additional O(i 2 ) terms. Not only that, knowing exactly how these v i, j coefficients look could aid in future questions asked by researchers interested in the M/M/c 2-class priority queue.
Readers should keep in mind that the expressions we have derived for the v i, j coefficients do contain binomial coefficients, and one should be careful to avoid roundoff errors while computing these expressions.

Theorem 3 (Coefficients)
The coefficients {v i, j } i≥ j≥0 from Theorem 2 are as follows: Our aim is to show v i+1, j also satisfies (4.1) for 0 ≤ j ≤ i + 1: We do this by substituting (4.1) into (3.12b) and simplifying. There are three cases to consider: (i) j = 1; (ii) 2 ≤ j ≤ i − 1; and (iii) j = i. We focus on case (ii), with cases (i) and (iii) following similarly.
We first examine the V 1 v i, j−1 term in (3.12b) by substituting (4.1). Here, Next, write (4.2) in a form where the binomial coefficients match the ones in (4.1): The remaining terms on the right-hand side of (3.12b) can be further simplified by substituting (4.1). Doing so reveals that Swapping the order of the triple summation in (4.4) gives The inner-most summation over k of (4.5) can be evaluated using Lemma 3 of Appendix A: Next, substitute (4.6) back into (4.5) and focus on the inner-most double summation of (4.5). This gives Substituting (4.7) into (4.5) and changing the two outer summation indices shows (4.8) Furthermore, since we can merge the single summation with the double summation in (4.8). In other words, (4.10) Finally, summing (4.3) and (4.10) produces (4.1), as Summing the coefficients in front of the binomial coefficient terms proves case (ii). Cases (i) and (iii) follow similarly.
We now have an explicit expression for the coefficients. Substitute the expressions for {v i, j } i≥ j≥0 into (3.11) to obtain, for j ≥ 1, Swapping the order of the double summation and grouping coefficients in front of each Laplace transform reveals the dependence of π (i,c−1+ j) (α) on the Laplace transforms π (0,c−1) (α), π (1,c−1) (α), . . . , π (i,c−1) (α): From this expression, we see that for each fixed i ≥ 0, as j → ∞, π (i,c−1+ j) (α) behaves in a manner analogous to that found in Theorem 3.1 of [23], which addresses, when c = 1, the asymptotic behavior of the stationary distribution as the number of high-priority customers approaches infinity, while the number of low-priority customers is fixed. The explicit expression (4.13) can be used to obtain an expression for the Laplace transforms of the number of class-1 customers in the system. That is, where we can simplify the final infinite sum as via the identity (4.16)

Laplace transforms for states in the horizontal boundary
In the previous section, we showed how to express each Laplace transform for the states on the vertical boundary and within the interior explicitly in terms of transforms for the states in the horizontal boundary. So, it remains to determine the transforms for the states in the horizontal boundary. In this section, we show that the latter Laplace transforms satisfy a variant of Ramaswami's formula, which will allow us to numerically compute these transforms recursively. The approach we use to compute the above-mentioned variant of Ramaswami's formula makes, like the CAP method, repeated use of Theorem 1. This approach is highly analogous to the approach used in [17] to study block-structured Markov processes, yet slightly modified since we are interested in recursively computing Laplace transforms only associated with states within the horizonal boundary. This idea of restricting ourselves to a subset of the state space seems similar in spirit to the censoring approach featured in the work of Li and Zhao [22], but it is not currently obvious to the authors if their approach is applicable to our setting.
We first introduce some relevant notation. Define the 1 × c row vectors π i (α) as To properly state the Ramaswami-like formula satisfied by these row vectors, we need to define additional matrices. First, we define the c × c transition rate submatrices corresponding to lower levels L i , i ≥ c, as A 1 := λ 1 I, A −1 := diag(cμ 1 , (c − 1)μ 1 , . . . , μ 1 ) and where I is the c × c identity matrix and diag(x) is a square matrix with the vector x along its main diagonal. We further define the c × c level-dependent transition rate submatrices associated with L i , 1 ≤ i ≤ c − 1, as A (i) Next, we define the collection of c × c matrices {W m (α)} m≥0 . Each element of W m (α) is equal to 0 except for element W m (α) c−1,c−1 , which is defined as We also need the collection of c × c matrices {G i, j (α)} i> j≥0 , where the (k, l)-th element of G i, j (α) is defined as Finally, we will need the collection of c × c matrices {N i (α)} i≥1 , whose elements are defined as follows: Note that N i (α) = N c (α) for i ≥ c, and we therefore denote N(α) := N c (α).

A Ramaswami-like recursion
The following theorem shows that the vectors of transforms {π i (α)} i≥0 satisfy a recursion analogous to Ramaswami's formula [25].
Proof This result can be proven by making use of the approach found in [17]. Using Theorem 1 with A = L ≤i , we see that for i ≥ 0 and 0 ≤ j ≤ c − 1, (5.6) Due to the structure of the transition rates, many terms in the summation of (5.6) are zero. In particular, (5.6) can be stated more explicitly as We now simplify each expectation appearing within the first sum on the right-hand side of (5.7). Summing over all ways in which the process reaches phase c − 1 again yields Applying the strong Markov property to each expectation appearing in (5.8) shows that where the last equality follows from the definitions of G l,i+1 (α) and N i+1 (α), and Lemma 8 of Appendix C. The expectations appearing within the second sum of (5.7) can easily be simplified by recognizing that they are elements of N i+1 (α). Hence, we ultimately obtain which, in matrix form, is (5.5).
It remains to derive computable representations of {G i, j (α)} i> j≥0 , as well as the matrices {N i (α)} 1≤i≤c .

Computing the G i, j (α) matrices
The next proposition shows that each G i,j (α) matrix can be expressed entirely in terms of the subset {G i+1,i (α)} 0≤i≤c−1 .

Proposition 1 For each pair of integers i, j satisfying i
Furthermore, for each integer k ≥ 0 we also have Proof Equation (5.11) can be derived by applying the strong Markov property in an iterative manner, while (5.12) follows from the homogeneous structure of X along all lower levels L i , i ≥ c.
In light of Proposition 1, our goal now is to determine the subset of matrices {G i+1,i (α)} 0≤i≤c−1 . We first focus on showing that G(α) := G c,c−1 (α) is the solution to a fixed-point equation. The next proposition shows that through successive substitutions one can obtain G(α) from (5.13). Now that we have a method for approximating G(α), it remains to find a method for computing G i+1,i (α), 0 ≤ i ≤ c − 2. The next proposition shows that these matrices can be computed recursively.

Proposition 4 For each integer i satisfying
Proof Similar to Proposition 2, (5.18) follows by a one-step analysis. The inverse in (5.18) exists because the matrix A (i+1) −1 is a diagonal matrix whose diagonal elements are all positive.

Computing the N i (α) matrices
The matrices {N i (α)} 1≤i≤c can be expressed in terms of {G i, j (α)} i≥ j≥0 .

Proposition 5 For each integer i satisfying
where we use the convention A (c) Proof The proof uses one-step analysis and is analogous to the proofs of Propositions 2 and 4. It is therefore omitted.

Computing π 0 (α)
It remains to devise a method for computing the vector π 0 (α) so that the Ramaswamilike recursion from Theorem 4 can be properly initialized. The following is an adaptation of [17,Section 3.3]. We define the c × c matrix N 0 (α) whose elements are given by In the derivation to follow, we require the notation (A) [i, j] which represents the matrix A with row i and column j removed (meaning it is a (c − 1) × (c − 1) matrix), while keeping the indexing of entries exactly as in A. Similarly, (A) [i,·] has row i removed from A (meaning it is a (c − 1) × c matrix) and (A) [·, j] has column j removed from A (meaning it is a c × (c − 1) matrix).

Proposition 6 We have
Proof Similar to the proofs of Propositions 2, 4 and 5, we use a one-step analysis and the strong Markov property to prove the result.
Remark 1 (An alternative method for determining π 0 (α)) The Kolmogorov forward equations can also be used to derive π 0 (α). The transition functions are known to satisfy these equations, since sup x∈S q(x) < ∞. Taking the Laplace transform of the Kolmogorov forward equations for the states in L 0 yields, after using (3.5), where e i is a row vector with all elements equal to zero except for the i-th element which is unity. Finally, using Theorem 4 to express π 1 (α) in terms of π 0 (α) yields If one is instead interested in the stationary distribution and in particular the stationary probabilities of L 0 , i.e., lim α↓0 α π 0 (α), (5.24) results in a homogeneous system of equations, but it is not clear that this system still has a unique solution. On the other hand, the approach outlined earlier for determining π 0 (α) can straightforwardly be employed to obtain lim α↓0 α π 0 (α).

Numerical implementation
In order to compute the vectors {π i (α)} i≥0 , we first need to compute the matrices The first step is to compute G(α) = G c,c−1 (α). Proposition 3 shows that this matrix can be approximated by using the recursion (5.16). Using this recursion requires us to truncate the infinite sum appearing within the recursion. One way of applying this truncation is as follows: given a fixed tolerance , pick an integer κ large enough so that Once κ has been found, we can use the approximation since the modulus of each element of the matrix on the left-hand side of (5.26) can be shown to be within of what is being approximated. Here we used that the matrices W l (α) only have one element and the absolute value of each element of Z(n, α) (and G(α)) is less than or equal to 1. Hence, we propose using the recursion to approximate G(α). Notice that we can determine κ satisfying (5.25) by writing the left-hand side of (5.25) as An explicit expression for the infinite sum on the right-hand side of (5.28) can be derived with the help of Lemma 8 of Appendix C and the generating function of the Catalan numbers; the finite sum can be computed numerically.
Once G(α) has been found, we can use Proposition 4 to compute each G i+1,i (α) matrix, for 0 ≤ i ≤ c − 2. For this computation, we use the same truncation procedure as outlined above.
Since we cannot compute this infinite sum, we determine ψ (i, j) (α) for all (i, j) in a sufficiently large bounding box S k := {(i, j) ∈ S : 0 ≤ i ≤ k} for some k ≥ 0. Notice that (4.15) allows S k to be an infinitely large rectangle. The choice of k in S k clearly influences the quality of the approximation. A simple procedure to choose k is the following. Define For c = 1 we can normalize the solution as outlined above, or we can explicitly determine the value of π (0,0) (α); see the next section.

The single-server case
We now turn our attention to the case where c = 1, i.e., the case where the system consists of a single server. In this case, the analysis of the Laplace transforms π (i,0) (α), i ≥ 0, simplifies considerably. The expressions for the Laplace transforms for the states in the interior and on the vertical boundary are identical to the multi-server case.
In light of (6.1), to evaluate π (0,0) (α) the only thing that needs to be determined is the expectation E (0,0) [e −ατ (0,0) ]. This quantity is the Laplace-Stieltjes transform of the sum of two independent exponential random variables: One is the exponential random variable E λ having rate λ, the other is the busy period B of an M/G/1 queue having arrival rate λ and hyperexponential service times having cumulative distribution function F(·). More specifically, The Laplace-Stieltjes transform ϕ(α) of B is known to satisfy the Kendall functional equation .

Proposition 7 The scalar G(α) is a solution to
Proof From Proposition 2, we easily find that when c = 1, The infinite series appearing in (6.8) can be simplified using Lemma 9 of Appendix C; doing so yields (6.7).
Even though we cannot use (6.7) to write down an explicit expression for G(α), we can still use it to devise an iterative scheme for computing G(α). The next result shows that N (α) can be expressed in terms of G(α).

Proposition 8 We have
. (6.9) Proof Using Proposition 5, we observe that when c = 1, The proof is completed by applying Lemma 9 of Appendix C to (6.10).
We now focus on the recursion for the horizontal boundary. When c = 1, Theorem 4 reduces to For the inner-most sum over l we have (6.12) Clearly, i + 1 − k ≥ 1. So let us try to evaluate the tail of the generating function of {w (6.13) which follows from an application of Lemma 9 of Appendix C. The remaining finite summation is easy to compute since each w (λ 2 ,μ 2 ,λ 1 ) m (α) term, by Lemma 8 of Appendix C, can be stated in terms of b K (·) functions, and these satisfy the recursion found in Lemma 4 of Appendix A. These observations allow us to state the following theorem, which yields a practical method for recursively computing Laplace transforms of the form π (i,0) (α).
Theorem 5 (Horizontal boundary, single server) When c = 1, the Laplace transforms of the transition functions on the horizontal boundary satisfy the following recursion: for i ≥ 0, (6.14)

Conclusion
In this paper, we analyzed an M/M/c priority system with two customer classes, class-dependent service rates and a preemptive resume priority rule. This queueing system can be modeled as a two-dimensional Markov process for which we analyzed the time-dependent behavior. More precisely, we obtained expressions for the Laplace transforms of the transition functions under the condition that the system is initially empty.
Using a slight modification of the CAP method, we showed that the Laplace transforms for the states with at least c high-priority customers can be expressed in terms of a finite sum of the Laplace transforms for the states with exactly c − 1 high-priority customers. This expression contained coefficients that satisfy a recursion. We solved this recursion to obtain an explicit expression for each coefficient. In doing so, each Laplace transform for the states on the vertical boundary and in the interior can easily be calculated from the values of the Laplace transforms for the states in the horizontal boundary.
Next, we developed a Ramaswami-like recursion for the Laplace transforms for the states in the horizontal boundary. The recursion required the collections of matrices {G i+1,i } 0≤i≤c−1 and {N i (α)} 1≤i≤c for which we showed that they can be determined iteratively. We demonstrated two ways in which the initial value of the recursion, i.e., the vector π 0 (α), can be calculated. Finally, we discussed the numerical implementation of our approach for the horizontal boundary.
In the single-server case, the expressions for the Laplace transforms for the states on the vertical boundary and in the interior were identical to the multi-server case. The expressions for the horizontal boundary, however, simplified considerably. Specifically, the initial value π (0,0) (α) of the recursion could be determined by comparing the queueing system to an M/G/1 queue with hyperexponentially distributed service times. Moreover, the calculation of G(α) and N (α), which are now scalars, simplified greatly.
We now comment on how our expressions for the Laplace transforms of the transitions functions can be used to determine the stationary distribution. It is clear from the transition rate diagram in Fig. 3 that the Markov process X is irreducible. Moreover, it is well-known that X is positive recurrent if and only if ρ < 1. In that case, X has a unique stationary distribution p := [p x ] x∈S . To compute each p x term from π x (α), simply note that p x = lim α↓0 α π x (α). (7.1) Using this observation, we see that the procedure for finding p is highly analogous to the one we presented for finding the Laplace transforms of the transition functions.
The rest of the proof follows from a direct application of [11,Chapter 2,Eq. (12.16)].

Lemma 3 We have the identity
Proof Change the summation variable to n = l − k to get proving the claim.
Proof The terms b K (z) appear in [5,Section 3.3], where the authors derive (A.9). We believe that a small typographical error appears in their recursion that we have fixed here. To do so, in [5], substitute (23) into (24) to obtain the correct form of (16). We now derive the recursion (A.10). Since we already have the explicit expression (A.8) for b K (z), we will substitute this into (A.10) and show that b K +1 (z) again is given by (A.8).
Rewrite the first term on the right-hand side of (A.10) as Substituting (A.8) into the finite sum of (A.10) gives Switch the order of the triple summation to obtain where we introduced n = l − k. Employing Lemma 2 for the inner-most summation results in (A.14) The double summation sums over all k, m ≥ 0 such that 0 ≤ k + m ≤ K . An equivalent summation is over the diagonals k + m = d with 0 ≤ d ≤ K and k, m ≥ 0: Using the identity C d+1 = d k=0 C k C d−k , setting k = d + 1 and rewriting the binomial coefficient produces Finally, summing (A.11) and (A.16) yields proving the recursion (A.10) is correct.

Appendix B: Random-product representation
The results we present in Appendix C make use of the random-product representation theory recently developed and discussed in [6,12]. Using this theory to study a given Markov process X requires selecting an additional Markov processX := {X (t)} t≥0 that shares the same state space S as X . The elements of the transition rate matrixQ ofX must satisfy the following two properties; see [6, Section 1] and [12, Section 2]: Associate with the Markov process X its transition times {T n } n≥0 , where T 0 := 0 and T n represents the n-th transition time of X . From the transition times, we create the embedded discrete-time Markov chain {X n } n≥0 as X n := X (T n ), n ≥ 0. The sequences {T n } n≥0 and {X n } n≥0 are constructed and defined similarly.
We will also need to make use of discrete-time hitting-time random variables. We define, for each set A ⊂ S, The random-product representation can be used to determine the Laplace transforms of the transition functions, when the process is assumed to start in a fixed state x ∈ S. We will employ the following theorem, which originally appeared in [ where .

Appendix C: Results for M/M/1 queues
Here we derive key quantities associated with M/M/1 queues.
Proof The time τ 0 is the minimum of the busy period of a regular M/M/1 queue, i.e., B λ,μ , and an exponential random variable with parameter θ . Note that the two random variables are independent. So, where E θ is an exponential random variable with parameter θ . Under this description, the time τ j is equal to the time τ * j to reach state j in a regular M/M/1 queue. Furthermore, conditioning on both τ * j and B λ,μ yields The remainder of the proof follows from the proof of [10, Lemma 4].
Now, change the summation index by setting l = k − 1 to retrieve so that Lemma 4 proves the claim (C.15).

Lemma 9
The generating function of the {w Proof We use the definition of w (λ,μ,θ) i (α) in Lemma 7 and condition on the length of the busy period: where we used the probability generating function of a Poisson distribution with parameter θ t.