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

Queueing Systems, Jul 2017

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs11134-017-9541-2.pdf

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

Time-dependent analysis of an M/ M/c preemptive priority system with two priority classes Jori Selen 0 1 2 Brian Fralix 0 1 2 B Jori Selen 0 1 2 0 Department of Mathematical Sciences, Clemson University , Clemson, SC , USA 1 Department of Mechanical Engineering, Eindhoven University of Technology , Eindhoven , Netherlands 2 Department of Mathematics and Computer Science, Eindhoven University of Technology , Eindhoven , Netherlands 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. Mathematics Subject Classification 60J27 · 60J35 · 60K25 Static priority; Time-dependent analysis; Laplace transforms; Multi-dimensional Markov process 1 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 are ρ1 = 1/3, ρ2 = 1/2 and c varies = 10−8. Parameter settings 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. ∞ ∞ 100 21 62 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. 2 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 class2 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 Xn(t ) represent the number of class-n customers in the system at time t , and define X (t ) := (X1(t ), X2(t )). Then, X := {X (t )}t≥0 is a CTMC on the state space S = N20. 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 q((i, j ), (i + 1, j )) = λ1, i, j ≥ 0, q((i, j ), (i, j + 1)) = λ2, i, j ≥ 0, q((i, j ), (i − 1, j )) = max(min(i, c − j ), 0)μ1, i ≥ 1, j ≥ 0, q((i, j ), (i, j − 1)) = min(c, j )μ2, i ≥ 0, j ≥ 1. X2(t) λ2 j2 cμ2 X1(t) Each transition function px,y (·) has a Laplace transform πx,y (·) that is well-defined on the subset of complex numbers C+ := {α ∈ C : Re(α) > 0} as 0 ∞ πx,y (α) := e−αt px,y (t ) dt, α ∈ C+. 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., px (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 . (2.1) (2.2) X2(t) vertical boundary c c − 1 2.1 Notation and terminology interior horizontal boundary c − 1 It helps to decompose the state space S into a countable number of levels, where, for each integer i ≥ 0, the i -th level is the set {(i, 0), (i, 1), . . .}. We further decompose the i -th level into an upper level and a lower level: Upper level i is defined as Ui := {(i, c), (i, c + 1), . . .}, while lower level i is simply Li := {(i, 0), (i, 1), . . . , (i, c − 1)} and the union of lower levels L0, L1, . . . , Li is denoted by L≤i = ik=0 Lk . The set of all states in phase j is denoted by Pj := {(0, j ), (1, j ), . . .}. We sometimes refer to upper level U0 as the vertical boundary. The union of upper levels is called the interior of the state space. Finally, the union U1 ∪ U2 ∪ · · · L0 ∪ L1 ∪ · · · 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 Ez [ f (Z )] represent the expectation of a functional of Z , conditional on Z (0) = z, and Pz (·) denotes the conditional probability associated with Ez [·]. In our analysis, it should be clear from the context what is being conditioned on when we write Pz (·) or Ez [·]. X1(t) (2.3) (2.4) 2.2 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 Pn(·), which, in this case, represents conditioning on Qλ,μ(0) = n, let Bλ,μ denote the busy period duration induced by these customers. Under P1(·), the Laplace– Stieltjes transform of Bλ,μ is given by φλ,μ(α) := E1[e−α Bλ,μ ] = λ + μ + α − (λ + μ + α)2 − 4λμ 2λ . Recall that under Pn(·), Bλ,μ is equal in distribution to the sum of n i.i.d. copies of Bλ,μ under the measure P1(·); see, for example, [28, p. 32]. Thus, for each integer n ≥ 1 we have En[e−α Bλ,μ ] = φλ,μ(α)n. 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 {Qλ,μ(t )}t≥0. For each integer i ≥ 0, define w(λ,μ,θ)(α) := E1[e−α Bλ,μ 1{Λθ (Bλ,μ) = i }]. i Lemma 7 of Appendix C shows that the w(λ,μ,θ)(α) terms satisfy a recursion, and, in i 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 We will also need to make use of hitting-time random variables. We define for each set A ⊂ S, τA := inf{t > 0 : lsi↑mt X (s) = X (t ) ∈ A} 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}. (2.5) and φ2 := φλ2,cμ2 (λ1 + α), r2 := ρ2φλ2,cμ2 (λ1 + α), ρ2φλ2,cμ2 (λ1 + α) Ω2 := λ2(1 − ρ2φλ2,cμ2 (λ1 + α)2) . Further results for M/M/1 queues are presented in Appendix C. (2.6) (2.7) (2.8) (2.9) (2.10) 2.3 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 Ui , 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 {vk,l }i≥k≥l≥0 that satisfy a recursion. 2. In Sect. 4, we obtain an explicit expression for the coefficients {vk,l }k≥l≥0. This in turn shows that for each i ≥ 0, each Laplace transform for each state in Ui , i.e., π(i,c−1+ j)(α), j ≥ 1, can be explicitly expressed in terms of π(k,c−1)(α), 0 ≤ k ≤ i . 3. 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 π(i, j)(α), i ≥ 0, 0 ≤ j ≤ c − 1. 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. 3 A slight modification of the CAP method The following theorem is used in multiple ways throughout our analysis. It appears in [17, Theorem 2.1] and can be derived by taking the Laplace transform of both sides of the equation at the top of page 124 of [ 21 ]. Equation (3.2) is the Laplace transform version of [10, Theorem 1]. Theorem 1 Suppose A and B are disjoint subsets of S with x ∈ A. Then for each y ∈ B, πx,y (α) = πx,z (α)(q(z) + α)Ez e−αt 1{X (t ) = y} dt , or, equivalently, z∈A πx,y (α) = πx,z (α) 3.1 Laplace transforms for states along the vertical boundary In this subsection, we employ Theorem 1 to express each Laplace transform π(0,c−1+ j)(α), j ≥ 1, in terms of π(0,c−1)(α). 0 τA 0 τA (3.1) (3.2) Using Theorem 1 with A = U0c 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 is equal in distribution to the minimum of the 0 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 can 0 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 λ2E(0,c) 0 τUc j 0 e−αt 1{X (t ) = (0, c − 1 + j )} dt = r2 . Substituting (3.4) into (3.3) then yields π(0,c−1+ j)(α) = π(0,c−1)(α)r2j , j ≥ 1. 3.2 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 Ui in terms of the transforms in upper level Ui−1 and in state (i, c − 1). Second, we use this result to express the transforms in upper level Ui in terms of the transforms for the states (0, c−1), (1, c−1), . . . , (i, c−1) and some additional coefficients. Employing again Theorem 1, now with A = Uic, yields, for i, j ≥ 1, 0 τUc i e−αt 1{X (t ) = (i, c − 1 + j )} dt 0 τUc i e−αt 1{X (t ) = (i, c − 1 + j )} dt π(i,c−1+ j)(α) = = z∈Uic ∞ k=1 τUc i e−αt 1{X (t ) = (i, c − 1 + j )} dt . (3.3) (3.4) (3.5) (3.6) The expectation λ1E(i,c−1+k) 0 where, for j, k ≥ 1, (3.11) (3.12a) (3.12b) 0 k=1(·), rep(3.9) Substituting (3.8) into (3.6) and simplifying yields a recursion. Specifically, for i ≥ 0 and j ≥ 1, Theorem 2 (Interior) For i ≥ 0, j ≥ 1, vi,k (1 − V2)k j − k1 + k r2j , vi+1,0 = π(i+1,c−1)(α), vi+1, j = V1 vi, j−1 + V2 vi,k , 1 ≤ j ≤ i + 1, i k= j 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 τUc i e−αt 1{X (t ) = (i, c − 1 + j )} dt = Υ ( j, k), j, k ≥ 1, (3.8) where the quantities {vi, j }i≥ j≥0 satisfy the following recursive scheme: for i ≥ 0, with initial condition v0,0 = π(0,c−1)(α). Here V1 = 1λ−1rΩ2φ22 , and V2 = r2φ2. Throughout we follow the convention that all empty sums, such as resent the number zero. Proof Clearly, when i = 0, (3.11) agrees with (3.5). Proceeding by induction, assume (3.11) holds among upper levels U0, U1, . . . , Ui for some i ≥ 0. Substituting (3.11) into (3.10) yields j − l1++1l + 1 r2j To further simplify the right-hand side of (3.14), increase the summation index of the first summation by one by setting k = l + 1, multiply its summands by 11−−rr22φφ22 , and change the order of the double summation. This yields j (3.14) = r2 π(i+1,c−1)(α) + V1vi,k−1(1 − V2)k j − k1 + k r2j + + i i V1V2 k=1 l=k j = r2 π(i+1,c−1)(α) i+1 k=1 V1 vi,k−1 + V2 vi,l (1 − V2)k j − k1 + k r2j i l=k vi,l (1 − V2)k j − k1 + k r2j , (3.15) which shows π(i+1,c−1+ j)(α) satisfies (3.11), completing the induction step. 4 Deriving an explicit expression for vi, j {vk,l }i≥k≥l≥0 a. Compute π(i+1,c−1)(α). b. Compute {vi+1,l }i+1≥l≥0. Theorem 2 suggests that the Laplace transforms for the states within the interior can be computed recursively: 1. Initialization step Determine π(0,c−1)(α), which yields each Laplace transform π(0,c−1+ j)(α), j ≥ 1. 2. Recursive step on i Given π(k,c−1)(α) for 0 ≤ k ≤ i and the coefficients 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 {vi,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 {vk,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 vi, 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 vi, j coefficients do contain binomial coefficients, and one should be careful to avoid roundoff errors while computing these expressions. Theorem 3 (Coefficients) The coefficients {vi, j }i≥ j≥0 from Theorem 2 are as follows: for i ≥ j ≥ 0, j vi, j = V1 π(i− j,c−1)(α) + i Proof From (3.12) we find, for each i ≥ 0, that vi,0 = π(i,c−1)(α) and vi,i = V1i π(0,c−1)(α); these expressions agree with (4.1). Next, assume for some integer i ≥ 0 that vi, j satisfies (4.1) for 0 ≤ j ≤ i . Our aim is to show vi+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 V1vi, j−1 term in (3.12b) by substituting (4.1). Here, j V1vi, j−1 = V1 π(i+1− j,c−1)(α) i + + i+1 Next, write (4.2) in a form where the binomial coefficients match the ones in (4.1): j (4.2) = V1 π(i+1− j,c−1)(α) + i+1 The remaining terms on the right-hand side of (3.12b) can be further simplified by substituting (4.1). Doing so reveals that V1V2 vi,k = k= j Swapping the order of the triple summation in (4.4) gives (4.4) = V k+1π(i−k,c−1)(α)V2 1 V l+1π(i−l,c−1)(α) 1 l− j l−m k m=1 k= j l − k The inner-most summation over k of (4.5) can be evaluated using Lemma 3 of Appendix A: l−m Next, substitute (4.6) back into (4.5) and focus on the inner-most double summation of (4.5). This gives i k= j l− j l−m k m=1 k= j l − k (4.6) (4.7) i+1 Substituting (4.7) into (4.5) and changing the two outer summation indices shows l V1 π(i+1−l,c−1)(α) m=2 l− j l + j m − j − m l − j m (l − j )(l − 1) we can merge the single summation with the double summation in (4.8). In other words, i+1 Finally, summing (4.3) and (4.10) produces (4.1), as V1vi, j−1 + V1V2 vi,k j = V1 π(i+1− j,c−1)(α) i k= j (4.8) (4.9) V2m . (4.10) (4.11) k= j+1 k= j+1 Summing the coefficients in front of the binomial coefficient terms proves case (ii). Cases (i) and (iii) follow similarly. V2m(1 − V2)k j − k1 + k r2j . (4.12) 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)(α): j π(i,c−1+ j)(α) = r2 π(i,c−1)(α) + V1l π(i−l,c−1)(α) (1 − V2)l j − l1 + l r2j 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, e−αt P(X1(t) = i | X (0) = (0, 0)) dt ∞ e−αt ∞ j=0 P(X (t) = (i, j ) | X (0) = (0, 0)) dt π(i, j)(α) = π(i, j)(α) + π(i,c−1+ j)(α), (4.14) 0 = = ∞ 0 ∞ j=0 where we can simplify the final infinite sum as ∞ j=1 + r2 π(i,c−1+ j)(α) = 1 − r2 π(i,c−1)(α) + l−1 (1 − V2)kr2 l−k k=0 (1 − r2)k+1 m=1 l − k k l=1 via the identity ∞ j=1 . (4.16) 5 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 π i (α) := π(i,0)(α) π(i,1)(α) · · · π(i,c−1)(α) , i ≥ 0. (5.1) 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 Li , i ≥ c, as A1 := λ1I, A−1 := diag(cμ1, (c − 1)μ1, . . . , μ1) and A0 := ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ −λ2 λ2 μ2 −(λ2 + μ2) 2μ2 λ2 −(λ2 + 2μ2) ⎤ ⎥ λ2 ⎥⎥ − A1 − A−1, . . . ⎥⎥ ⎦ (c − 1)μ2 −(λ2 + (c − 1)μ2) have A(0) 0 := A0 + A−1. 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 Li , 1 ≤ i ≤ c − 1, as A(i) −1 := diag(x(i)) with (x(i)) j := min(i, c − j )μ1, 0 ≤ j ≤ c − 1, A(i) 0 := A0 + A−1 − A(−i)1, and for L0 we (5.2) Next, we define the collection of c × c matrices {Wm (α)}m≥0. Each element of Wm (α) is equal to 0 except for element Wm (α) c−1,c−1, which is defined as Gi, j (α) k,l := E(i,k)[e−ατL j 1{X (τL j ) = ( j, l)}], 0 ≤ k, l ≤ c − 1. Finally, we will need the collection of c × c matrices {Ni (α)}i≥1, whose elements are defined as follows: Ni (α) k,l := E(i,k) 0 τLi−1 e−αt 1{X (t ) = (i, l)} dt , 0 ≤ k, l ≤ c − 1. Note that Ni (α) = Nc(α) for i ≥ c, and we therefore denote N(α) := Nc(α). 5.1 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 ]. Theorem 4 (Horizontal boundary) For each integer i ≥ 0, we have where we use the convention Gi+1,i+1(α) = I. 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, l=i+1 Wl−k (α)Gl,i+1(α)Ni+1(α), (5.5) (5.3) (5.4) (5.7) E(k,c) = E(k,c) 0 = ∞ l=i+1 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 = ∞ l=i+1 ∞ l=i+1 i + m=0 (5.8) (5.9) where the last equality follows from the definitions of Gl,i+1(α) and Ni+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 Ni+1(α). Hence, we ultimately obtain l=i+1 π(i,m)(α) A1 m,m Ni+1(α) m, j , (5.10) which, in matrix form, is (5.5). It remains to derive computable representations of {Gi, j (α)}i> j≥0, as well as the matrices {Ni (α)}1≤i≤c. 5.2 Computing the Gi, j (α) matrices of the subset {Gi+1,i (α)}0≤i≤c−1. The next proposition shows that each Gi,j(α) matrix can be expressed entirely in terms Proposition 1 For each pair of integers i, j satisfying i > j ≥ 0, we have Gi, j (α) = Gi,i−1(α)Gi−1,i−2(α) · · · G j+1, j (α). Furthermore, for each integer k ≥ 0 we also have Gc+k,c−1+k (α) = Gc,c−1(α). 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 Li , i ≥ c. In light of Proposition 1, our goal now is to determine the subset of matrices {Gi+1,i (α)}0≤i≤c−1. We first focus on showing that G(α) := Gc,c−1(α) is the solution to a fixed-point equation. Proposition 2 The matrix G(α) satisfies G(α) = αI − A0 − W0(α) −1 A−1 + A1G(α)2 + Wl (α)G(α)l+1 . (5.13) 0 τ(Lc∪Uc)c Proof Assuming the matrix αI − A0 − W0(α) is invertible, Eq. (5.13) follows from a one-step analysis argument. To show that αI − A0 − W0(α) is indeed invertible, define the c × c matrix H(α) with elements H(α) i, j := E(c,i) e−αt 1{X (t ) = (c, j )} dt , 0 ≤ i, j ≤ c − 1, (5.14) and use again one-step analysis to establish ∞ l=1 αI − A0 − W0(α) H(α) = I, which proves the claim. The next proposition shows that through successive substitutions one can obtain G(α) from (5.13). Proposition 3 Suppose the sequence of matrices {Z(n, α)}n≥0 satisfies the recursion Z(n + 1, α) = αI − A0 − W0(α) −1 with initial condition Z(0, α) = 0. Then, (5.11) (5.12) (5.15) (5.16) (5.17) Proof This proof makes use of Proposition 2, and is completely analogous to the proofs of [18, Theorems 3.1 and 4.1] and [17, Theorem 3.4]. It is therefore omitted. Now that we have a method for approximating G(α), it remains to find a method for computing Gi+1,i (α), 0 ≤ i ≤ c − 2. The next proposition shows that these matrices can be computed recursively. Proposition 4 For each integer i satisfying 0 ≤ i ≤ c − 2, we have Gi+1,i (α) = αI − A(0i+1) − A1Gi+2,i+1(α) − ∞ l=i+1 Wl−(i+1)(α)Gl,i+1(α) −1A(−i +11). (5.18) Proof Similar to Proposition 2, (5.18) follows by a one-step analysis. The inverse in (5.18) exists because the matrix A(−i +11) is a diagonal matrix whose diagonal elements are all positive. We now have an iterative procedure for computing all {Gi+1,i (α)}0≤i≤c−1 matrices: first compute G(α) from Proposition 3, then use Proposition 4 to compute Gc−1,c−2(α), then Gc−2,c−3(α), and so on, stopping at G1,0(α). 5.3 Computing the Ni(α) matrices The matrices {Ni (α)}1≤i≤c can be expressed in terms of {Gi, j (α)}i≥ j≥0. Proposition 5 For each integer i satisfying 1 ≤ i ≤ c, we have Ni (α) = αI − A(i) 0 − A1Gi+1,i (α) − Wl−i (α)Gl,i (α) −1 , (5.19) ∞ l=i where we use the convention A(c) 0 = A0. Proof The proof uses one-step analysis and is analogous to the proofs of Propositions 2 and 4. It is therefore omitted. 5.4 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 N0(α) whose elements are given by N0(α) i, j := E(0,i) e−αt 1{X (t ) = (0, j )} dt , 0 ≤ i, j ≤ c − 1. (5.20) 0 τ(0,0) 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. We employ (3.2) and Proposition 5 to determine the elements of the row vector π 0(α). Set A = {(0, 0)} in Theorem 1: then, for 1 ≤ j ≤ c − 1 and c ≥ 2, π(0, j)(α) = = π(0,0)(α) λ1 G1,0(α) 0,l N0(α) l, j + λ2 N0(α) 1, j . (5.22) Hence, the transforms π(0, j)(α), 1 ≤ j ≤ c−1, can be expressed in terms of π(0,0)(α). In Sect. 5.5, we describe a numerical procedure to determine π(0,0)(α). 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 supx∈S q(x ) < ∞. Taking the Laplace transform of the Kolmogorov forward equations for the states in L0 yields, after using (3.5), π 0(α) αI − A(00) − π 0(α)W0(α) − π 1(α)A(−11) = e0, where ei 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 −π 0(α) −αI + A(0) 0 + W0(α) + A1 + Wl (α)Gl,1(α) N1(α)A(−11) = e0. (5.23) (5.24) If one is instead interested in the stationary distribution and in particular the stationary probabilities of L0, 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(α). 5.5 Numerical implementation In order to compute the vectors {π i (α)}i≥0, we first need to compute the matrices {Gi+1,i (α)}0≤i≤c−1, {Ni (α)}1≤i≤c, and N0(α) [0,0]. The first step is to compute G(α) = Gc,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 ∞ l=κ +1 |wl(λ2,cμ2,λ1)(α)| ≤ /λ2. Once κ has been found, we can use the approximation κ l=1 Wl (α)Z(n, α)l+1 ≈ Wl (α)Z(n, α)l+1, 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 Wl (α) 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 Z(n + 1, α) = αI − A0 − W0(α) −1 l=1 κ l=1 κ l=1 (5.25) (5.26) to approximate G(α). Notice that we can determine κ satisfying (5.25) by writing the left-hand side of (5.25) as ∞ l=κ +1 ∞ l=1 |wl(λ2,cμ2,λ1)(α)| = |wl(λ2,cμ2,λ1)(α)| − |wl(λ2,cμ2,λ1)(α)|. (5.28) 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. where π(0,0)(α) is an unknown transform and ψ(0,0)(α) = 1. It is clear that the ψx (α) can be computed using the same procedure as for πx (α), and (5.22) shows that ψ(0,0)(α), ψ(0,1)(α), . . . , ψ(0,c−1)(α) are computable expressions. We can calculate π(0,0)(α) from the normalization condition which yields x∈S πx (α) = π(0,0)(α) π(0,0)(α) = α 1 x∈S ψx (α) . Since we cannot compute this infinite sum, we determine ψ(i, j)(α) for all (i, j ) in a sufficiently large bounding box Sk := {(i, j ) ∈ S : 0 ≤ i ≤ k} for some k ≥ 0. Notice that (4.15) allows Sk to be an infinitely large rectangle. The choice of k in Sk clearly influences the quality of the approximation. A simple procedure to choose k is the following. Define k := x∈Sk ψx (α). Pick small and positive and continue increasing k until | k+|1−k| k | < . Then, set π(0,0)(α) = 1/(α k+1) to normalize the Laplace transforms πx (α) = π(0,0)(α)ψx (α). 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. 6 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. (5.29) (5.30) (5.31) (5.32) From [12, Corollary 2.1]—see also Theorem 6 in Appendix B—we have π(0,0)(α) = (q((0, 0)) + α)(1 − E(0,0)[e−ατ(0,0) ]) . 1 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, F (t ) = λ1 (1 − e−μ1t ) + λλ2 (1 − e−μ2t ), t ≥ 0. λ The Laplace–Stieltjes transform ϕ(α) of B is known to satisfy the Kendall functional equation . Furthermore, ϕ(α) can be determined numerically through successive substitutions of (6.3), starting with ϕ(α) = 0; see [1, Section 1] for details. Using independence of Eλ and B, E(0,0)[e−ατ(0,0) ] = E[e−α(Eλ+B) meaning that (see also [2, Eq. (36)]) for c = 1, λ ϕ(α), ] = λ + α 1 π(0,0)(α) = λ(1 − ϕ(α)) + α . We now turn our attention to the horizontal boundary. When c = 1, the matrices G(α) and N(α) become scalars, which we denote as G(α) and N (α), respectively. More precisely, G(α) := E(i+1,0)[e−ατLi ], N (α) := E(i,0) which are independent of i ≥ 1. Proposition 7 The scalar G(α) is a solution to (λ + μ1 + α)G(α) = μ1 + λ1G(α)2 + λ2G(α)φλ2,μ2 (λ1(1 − G(α)) + α). (6.7) Proof From Proposition 2, we easily find that when c = 1, (λ + μ1 + α)G(α) = μ1 + λ1G(α)2 + λ2 w(λ2,μ2,λ1)(α)G(α)l+1. l (6.8) τLi−1 e−αt 1{X (t ) = (i, 0)} dt , 0 ∞ l=0 (6.1) (6.2) (6.3) (6.4) (6.5) (6.6) 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 N (α) = α + λ + μ1 − λ1G(α) − λ2φλ2,μ2 (λ1(1 − G(α)) + α) . Proof Using Proposition 5, we observe that when c = 1, N (α) = α + λ + μ1 − λ1G(α) − λ2 wl(λ2,μ2,λ1)(α)G(α)l −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 π(i+1,0)(α) = λ1π(i,0)(α)N (α) i + λ2 ∞ For the inner-most sum over l we have w(λ2,μ2,λ1)(α)G(α)l−(i+1) N (α). l−k (6.11) 1 ∞ (6.9) (6.10) wl(−λ2k,μ2,λ1)(α)G(α)l−(i+1) = G(α)k−(i+1) wm(λ2,μ2,λ1)(α)G(α)m . (6.12) Clearly, i + 1 − k ≥ 1. So let us try to evaluate the tail of the generating function of {wm(λ2,μ2,λ1)(α)}m≥0 evaluated at the point G(α), i.e., ∞ m=i+1−k ∞ m=i+1−k wm(λ2,μ2,λ1)(α) G(α)m = φλ2,μ2 (λ1(1 − G(α)) + α) − wm(λ2,μ2,λ1)(α) G(α)m , (6.13) which follows from an application of Lemma 9 of Appendix C. The remaining finite summation is easy to compute since each wm(λ2,μ2,λ1)(α) term, by Lemma 8 of Appendix C, can be stated in terms of bK (·) 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, π(i+1,0)(α) = λ1π(i,0)(α)N (α) + λ2 i π(k,0)(α)G(α)k−(i+1) φλ2,μ2 (λ1(1 − G(α)) + α) w(λ2,μ2,λ1)(α)G(α)l N (α). l (6.14) 7 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 {Gi+1,i }0≤i≤c−1 and {Ni (α)}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 := [ px ]x∈S. To compute each px term from πx (α), simply note that px = lαi↓m0 α πx (α). 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. Acknowledgements The authors thank both the editor and the referees for their constructive comments on a previous version of our manuscript, as these helped us improve both the content and the readability of our paper. The first author is supported by a free competition Grant from the Netherlands Organisation for Scientific Research (NWO), while the second author gratefully acknowledges the support of the National Science Foundation (NSF), via Grant NSF-CMMI-1435261. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Appendix A: Combinatorial identities In this appendix, we collect some combinatorial identities that are used throughout the paper. These lemmas are likely known, but we prove them here to make the paper self-contained. Lemma 1 For j ≥ 1, l ≥ 0 and Υ ( j, k) defined in (3.9), we have (7.1) = λ1Ω2 k − 1 + l l r2k r j j − 1 + l + 1 l + 1 λ1Ω2 2 + 1 − r2φ2 r2φ2 1 (1 − r2φ2)l−m j − 1 + m m l Proof The result is nearly identical to [10, Lemma 1], so we omit its proof. Lemma 2 We have the identity K k=0 (K − k + 2l)! (k + 2m)! (K − k)! K + 2l + 2m + 1 2l + 2m + 1 . Proof Rewriting the fractions as two binomial coefficients, (A.2) = (2l)!(2m)! K k=0 K − k + 2l K − k k + 2m k . (A.1) (A.2) (A.3) The rest of the proof follows from a direct application of [11, Chapter 2, Eq. (12.16)]. Lemma 3 We have the identity l−m Proof Change the summation variable to n = l − k to get Splitting the summation and using the identity l−m proving the claim. Lemma 4 Define (A.4) (A.5) (A.6) (A.7) (A.8) (A.9) (A.10) n=m n m l l− j−1 = m n=m−1 l(m + 1) = m(l + 1 − j ) n m l− j l l− j = m n=m n − 1 m − 1 − n m l l − j = m m l + 1 − j , m + 1 l− j n m bK (z) := Ck K k=0 K + k zk , K − k 1 2k are the Catalan numbers. The sequence {bK (z)}K ≥0 satisfies where Ck := k+1 k two recursions. For K ≥ 2 it satisfies (K + 1)bK (z) = (2K − 1)(1 + 2z)bK −1(z) − (K − 2)bK −2(z), alternatively, for K ≥ 0, it satisfies bK +1(z) = bK (z) + z with b0(z) = 1 and b1(z) = 1 + z. Proof The terms bK (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 bK (z), we will substitute this into (A.10) and show that bK +1(z) again is given by (A.8). Rewrite the first term on the right-hand side of (A.10) as bK (z) = K k=0 Ck K + k K − k zk = K k=0 (A.11) Substituting (A.8) into the finite sum of (A.10) gives (A.13) (A.14) (A.15) Switch the order of the triple summation to obtain z K l=0 bl (z)bK −l (z) = z K K −k k=0 m=0 K K −k k=0 m=0 K l l=0 k=0 K l K −l l=0 k=0 m=0 Ck Cm zk+m+1 Ck Cm zk+m+1 Ck ll −+ kk zk K −l K − l + m zm K − l − m l + k Ck Cm l − k K − l + m K − l − m (A.12) K −m l=k K −k−m K − l + m K − l − m K − k − m − n + 2m , K − k − m − n where we introduced n = l − k. Employing Lemma 2 for the inner-most summation results in (A.13) = K K −k k=0 m=0 Ck Cm zk+m+1 K + 1 + k + m . K − (k + m) 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: Ck Cm zk+m+1 K + 1 + k + m K − (k + m) d=0 k,m≥0 : k+m=d K d = d=0 k=0 Ck Cd−k zd+1 K + 1 + d . K − d Using the identity Cd+1 = binomial coefficient produces (A.15) = K K d=0 K +1 Finally, summing (A.11) and (A.16) yields bK (z) + z 2k K + 1 + k Ck = K +1 Ck K + k K + 1 − k zk (A.16) (A.17) K k=0 K +1 k=0 K +1 Ck K + 1 + k K + 1 − k zk = bK +1(z), 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 process X˜ := {X˜ (t )}t≥0 that shares the same state space S as X . The elements of the transition rate matrix Q˜ of X˜ must satisfy the following two properties; see [6, Section 1] and [12, Section 2]: (i) For each x ∈ S, q˜ (x ) := −q˜ (x , x ) = y=x q˜ (x , y) = q(x ); (ii) For each x , y ∈ S, x = y, q˜ (x , y) > 0 if and only if q(y, x ) > 0. Associate with the Markov process X its transition times {Tn}n≥0, where T0 := 0 and Tn represents the n-th transition time of X . From the transition times, we create the embedded discrete-time Markov chain {Xn}n≥0 as Xn := X (Tn), 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, ηA := inf{n ≥ 1 : Xn ∈ A} as the first time the embedded chain make a transition into A. ηx should be understood to mean η{x}. The continuous-time hitting-time random variable τ˜A is defined anal(B.1) ogously to the definition in Sect. 2.1, and η˜ A represents the first time the embedded DTMC of X˜ reaches the set A. 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 [12, Corollary 2.1 and Theorem 2.1]. Theorem 6 Suppose y ∈ S, where y = x . Then the Laplace transform πx,y (α) of px,y (·) satisfies where πx,y (α) = πx,x (α)Ey e−ατ˜x η˜x q(X˜ l , X˜ l−1) , Appendix C: Results for M/M/1 queues Here we derive key quantities associated with M/M/1 queues. The first lemma is a restatement of [10, Lemma 4], which was inspired by Problems 22 and 23 of [19, Chapter 7]. Lemma 5 Suppose {X (t )}t≥0 is an M/M/1 queueing model with exponential clearings. Arrivals occur according to a Poisson process with rate λ, each service is exponentially distributed with rate μ and there is an external Poisson process having rate θ of clearing instants, where, whenever a clearing occurs, all customers in the system at the clearing time are removed. Then, for j ≥ 1, E1 1{τ j < τ0}e−ατ j = μλ φλ,μ(θ + α) j−1(1 − μλ φλ,μ(θ + α)2) 1 − μλ φλ,μ(θ + α)2 j . (C.1) 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, E1 1{τ j < τ0}e−ατ j = E1 1{τ j < Bλ,μ}1{τ j < Eθ }e−ατ j , (C.2) 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 (C.2) = E1 1{τ ∗j < Bλ,μ}1{τ ∗j < Eθ }e−ατ ∗j = E1 E 1{τ ∗j < Bλ,μ}1{τ ∗j < Eθ }e−ατ ∗j | Bλ,μ, τ ∗j = E1 1{τ ∗j < Bλ,μ}e−ατ ∗j E 1{τ ∗j < Eθ } | Bλ,μ, τ ∗j = E1 1{τ ∗j < Bλ,μ}e−(θ+α)τ ∗j . 0 0 τ0 0 τ0 The remainder of the proof follows from the proof of [10, Lemma 4]. The following lemma is a minor generalization of [10, Theorem 2], in that we verify it is still valid for α ∈ C+. Lemma 6 Suppose {X (t )}t≥0 is the clearing model of Lemma 5. Then for each j, k ≥ 1, Ek 0 τ0 ⎪⎨⎧ λ(1−μλ μφλλφ,λμ,(μθ(+θ+α)α)2) μλ φλ,μ(θ + α) j−k 1 − μλ φλ,μ(θ + α)2 k , 1 ≤ k ≤ j − 1, k ≥ j. Proof Define the Laplace transform of the transition functions π0, j (α) := e−αt p0, j (t ) dt. The state space for this single-server queue is S = N0. Select A = {0} and B = S\ A and apply Theorem 1 to obtain, for j ≥ 1, π0, j (α) = π0,0(α)(λ + α)E0 e−αt 1{X (t ) = j } dt . We can use the random-product representation of Appendix B to derive another expression for π0, j (α). Construct a Markov process X˜ := {X˜ (t )}t≥0 with transition rates q˜ (i, i − 1) = λ + θ and q˜ (i, i + 1) = μ for i ≥ 1. X˜ also has transitions from state 0 to every other state, but these do not factor into the calculations so there is no need to formally define them here. From Theorem 6, π0, j (α) = π0,0(α)E j e−ατ˜0 = π0,0(α)E1 e−α Bμ,λ+θ η˜0 q(X˜ l , X˜ l−1) l=1 q˜ (X˜ l−1, X˜ l ) λ λ + θ Dμ,λ+θ j (C.3) (C.4) (C.5) (C.6) λ μ = π0,0(α) φλ,μ(θ + α) j . Combining (C.6) with (C.7) gives E1 0 τ0 φλ,μ(θ + α) j . For now, abbreviate φ := φλ,μ(θ +α) and r := μλ φλ,μ(θ +α). The remaining expected values can be computed. First, for 2 ≤ k ≤ j , e−αt 1{X (t ) = j } dt = E1 1{τk < τ0} τ0 = E1 1{τk < τ0}e−ατk e−α(t−τ0)1{X (t ) = j } dt = E1 1{τk < τ0}e−ατk Ek τ0 τk 0 τ0 (C.8) (C.10) (C.11) (C.12) Furthermore, from Lemma 5, E1 0 τ0 meaning 0 τ0 E1 1{τk < τ0}e−ατk = r k−1(1 − r φ) 1 − (r φ)k , Ek 0 τ0 e−αt 1{X (t ) = j } dt = λ(1 −r r φ) r j−k (1 − (r φ)k ). Deriving the expected values when k > j is a little more straightforward. Here, Ek e−αt 1{X (t ) = j } dt = Ek 1{τ j < τ0}e−ατ j e−α(t−τ j )1{X (t ) = j } dt which proves the claim. sion Lemma 7 The expectation w(λ,μ,θ)(α), defined in (2.8), satisfies for i ≥ 1 the recuri θ λ w(λ,μ,θ)(α) = λ+μ+θ +α i 1 i w(−λ,μ,θ)(α)+ λ+μ+θ +α wi(λk,μ,θ)(α)wk(λ,μ,θ)(α) − (C.13) = Ek 1{τ j < τ0}e−ατ j E j e−αt 1{X (t ) = j } dt φk− j (1 − (r φ) j ), τ0 τ j 0 τ0 k=0 and w(λ,μ,θ)(α) = φλ,μ(θ + α). 0 Proof Conditioning on the length of the busy period, we have, for i = 0, w(λ,μ,θ)(α) = E1[e−α Bλ,μ 1{Λθ (Bλ,μ) = 0}] = E1[e−α Bλ,μ 1{Bλ,μ < Eθ }] 0 0 P(t < Eθ )e−αt f Bλ,μ (t ) dt = φλ,μ(θ + α). The recursion for i ≥ 1 follows from a one-step analysis and the strong Markov property. Since the birth–and–death process starts with 1 customer, the first event occurs after Eλ+μ+θ time and is either an arrival according to the Poisson process with probability θ /(λ + μ + θ ) or an arrival of an additional customer with probability λ/(λ + μ + θ ). If the former occurs, due to the strong Markov property, one less Poisson point needs to arrive. If the latter occurs, exactly i Poisson points need to arrive in two busy periods (due to the homogeneous structure of the birth–and–death process). This reasoning establishes the recursion (C.13). Lemma 8 The {wi(λ,μ,θ)(α)}i≥0 of Lemma 7 are given by w(λ,μ,θ)(α) = φλ,μ(θ + α) 0 and for i ≥ 1, i−1 k=0 w(λ,μ,θ)(α) = W1i φλ,μ(θ + α) i W k 2 = W1i φλ,μ(θ + α)bi−1(W2), W1 = λ(1 − 2φλ,μ(θ + α)) + μ + θ + α λφλ,μ(θ + α) W2 = λ(1 − 2φλ,μ(θ + α)) + μ + θ + α , . Proof For brevity, define wi := w(λ,μ,θ)(α). Rewrite (C.13) as i Using w0 = φλ,μ(θ + α), this reduces to (λ + μ + θ + α)wi+1 = θ wi + λ wi+1−k wk + 2w0wi+1 . (C.17) Straightforwardly substituting (C.15) into (C.18) and dividing by W i+1φλ,μ(θ + α) 1 results in bi (W2) = bi−1(W2) + W2 bi−k (W2) bk−1(W2). wi+1−k wk . (C.18) θ i k=1 i k=1 (C.16) (C.19) Now, change the summation index by setting l = k − 1 to retrieve bi−1−l (W2) bl (W2), l=0 so that Lemma 4 proves the claim (C.15). Lemma 9 The generating function of the {wi(λ,μ,θ )(α)}i≥0 is, for |z| < 1, i=0 wi(λ,μ,θ )(α) zi = φλ,μ(θ (1 − z) + α). Proof We use the definition of w(λ,μ,θ )(α) in Lemma 7 and condition on the length i of the busy period: E1 e−α Bλ,μ 1{Λθ ( Bλ,μ) = i } zi i=0 i=0 0 0 0 E1 e−α Bλ,μ 1{Λθ ( Bλ,μ) = i }zi = E1 e−α Bλ,μ zΛθ (Bλ,μ) E1 e−α Bλ,μ zΛθ (Bλ,μ) | Bλ,μ = t f Bλ,μ (t ) dt e−αt E zΛθ (t) f Bλ,μ (t ) dt e−(θ (1−z)+α)t f Bλ,μ (t ) dt = φλ,μ(θ (1 − z) + α), (C.22) where we used the probability generating function of a Poisson distribution with parameter θ t . 1. Abate, J., Whitt, W.: Solving probability transform functional equations for numerical inversion. Oper. Res. Lett. 12(5), 275–281 (1992) 2. Abate, J., Whitt, W.: Transient behavior of the M/G/1 workload process. Oper. Res. 42(4), 750–764 (1994) 3. Abate, J., Whitt, W.: Numerical inversion of Laplace transforms of probability distributions. ORSA J. Comput. 7(1), 36–43 (1995) 4. Abate, J., Whitt, W.: A unified framework for numerically inverting Laplace transforms. INFORMS J. Comput. 18(4), 408–421 (2006) 5. Abate , J. , Whitt , W. : Integer sequences from queueing theory . J. Integer Seq . 13 ( 5 ), 1 - 21 ( 2010 ) 6. Buckingham , P. , Fralix , B.: Some new insights into Kolmogorov's criterion, with applications to hysteretic queues . Markov Process. Relat. Fields 21 ( 2 ), 339 - 368 ( 2015 ) 7. Cobham , A. : Priority assignment in waiting line problems . Oper. Res . 2 ( 1 ), 70 - 76 ( 1954 ) 8. Davis , R.: Waiting-time distribution of a multi-server, priority queuing system . Oper. Res . 14 ( 1 ), 133 - 136 ( 1966 ) 9. den Iseger, P.: Numerical transform inversion using Gaussian quadrature . Prob. Eng. Inf. Sci. 20 , 1 - 44 ( 2006 ) 10. Doroudi , S. , Fralix , B. , Harchol-Balter , M. : Clearing analysis on phases: exact limiting probabilities for skip-free, unidirectional, quasi-birth-death processes ( 2015 ). ArXiv preprint arXiv:1503.05899v3 11. Feller , W.: An Introduction to Probability Theory and Its Applications, vol. I, 3 revised edn . Wiley, New York ( 1968 ) 12. Fralix , B. : When are two Markov chains similar? Stat . Prob. Lett. 107 , 199 - 203 ( 2015 ) 13. Gail , H. , Hantler , S. , Taylor , B.: On a preemptive Markovian queue with multiple servers and two priority classes . Math. Oper. Res . 17 ( 2 ), 365 - 391 ( 1992 ) 14. Harchol-Balter , M. , Osogami , T. , Scheller-Wolf , A. , Wierman , A. : Multi-server queueing systems with multiple priority classes . Queueing Syst . 51 ( 3-4 ), 331 - 360 ( 2005 ) 15. Jaiswal , N.: Preemptive resume priority queue . Oper. Res . 9 ( 5 ), 732 - 742 ( 1961 ) 16. Jaiswal , N.: Priority Queues . Academic Press Inc, New York ( 1968 ) 17. Joyner , J. , Fralix , B. : A new look at block-structured Markov processes . In: Working Paper ( 2016 ). http://bfralix.people.clemson.edu/preprints/BlockStructuredPaper8June.pdf 18. Joyner , J. , Fralix , B. : A new look at Markov processes of G/M/1-type . Stoch. Models 32 ( 2 ), 253 - 274 ( 2016 ) 19. Karlin , S. , Taylor , H.: A First Course in Stochastic Processes, 2nd edn . Academic Press, San Diego ( 1975 ) 20. Katehakis , M. , Smit , L. , Spieksma , F. : A comparative analysis of the successive lumping and the lattice path counting algorithms . J. Appl. Prob . 53 ( 1 ), 106 - 120 ( 2016 ) 21. Latouche , G. , Ramaswami , V. : Introduction to Matrix Analytic Methods in Stochastic Modeling . Society for Industrial and Applied Mathematics , Philadelphia ( 1999 ) 22. Li , Q.L. , Zhao , Y.Q. : The RG-factorization in block-structured Markov renewal processes . In: Zhu, X . (ed.) Observation, Theory and Modeling of Atmospheric Variability, pp. 545 - 568 . World Scientific, Singapore ( 2004 ) 23. Li , H. , Zhao , Y. : Exact tail asymptotics in a priority queue-characterizations of the preemptive model . Queueing Syst . 63 ( 1-4 ), 355 - 381 ( 2009 ) 24. Miller , D. : Computation of steady-state probabilities for M/M/1 priority queues . Oper. Res . 29 ( 5 ), 945 - 958 ( 1981 ) 25. Ramaswami , V.: A stable recursion for the steady state vector in Markov chains of M/G/1 type . Stoch. Models 4 ( 1 ), 183 - 188 ( 1988 ) 26. Sleptchenko , A., van Harten , A ., van der Heijden, M.: An exact solution for the state probabilities of the multi-class, multi-server queue with preemptive priorities . Queueing Syst . 50 ( 1 ), 81 - 107 ( 2005 ) 27. Sleptchenko , A. , Selen , J. , Adan , I., van Houtum, G. : Joint queue length distribution of multi-class, single-server queues with preemptive priorities . Queueing Syst . 81 ( 4 ), 379 - 395 ( 2015 ) 28. Takács , L. : Introduction to the Theory of Queues . Oxford University Press Inc, New York ( 1962 ) 29. Wang , J. , Baron , O. , Scheller-Wolf , A.: M/M/c queue with two priority classes . Oper. Res . 63 ( 3 ), 733 - 749 ( 2015 )


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs11134-017-9541-2.pdf

Jori Selen, Brian Fralix. Time-dependent analysis of an M / M / c preemptive priority system with two priority classes, Queueing Systems, 2017, 1-37, DOI: 10.1007/s11134-017-9541-2