Introduction

Matrix-analytic method (MAM) is an established powerful technique for the analysis of stochastic models which uses the special structure of a Markov chain to obtain its steady-state or transient distribution [15,11]. The basic step of a stochastic model analysis is to obtain the stability conditions of the model. The stability conditions within the MAM framework [15] are usually a version [20] of somewhat intuitive negative drift criterion [9,13]. However, the structure of a model may have additional properties that allow one to express the stability condition in an explicit form in terms of the model parameters, such as the input, service rates etc. Such properties are in the focus of this paper.

Consideration is given to multiclass multiserver queueing systems with the specific feature: automata-type transitions. Their structure allows to make a certain Markov chain embedding and, under certain assumptions, obtain the stability criterion involving only a product form (see (11) below). Surprisingly, it turns out that such a product form may be obtained even in the cases when some transition rates are not known (see (4) and discussion below).

This paper is inspired by the results on stability of a simultaneous service multiserver model describing a supercomputer, which were obtained in [17]. The model itself is of a separate interest. It is a multiserver model with twodimensional customers, namely, each such a customer requires a random number of servers for a (same at each occupied server) random time, while servers dedicated to a customer are seized and released simultaneously. Such a model, to

Stability Criterion of a QBD Process

The so-called quasi-birth-death process (QBD process) is a two-dimensional discrete-state continuous-time Markov process {X(t), Y (t)} t 0 , such that the level, X(t), is a nonnegative integer-valued variable increased or decreased by at most one at each transition, whereas the phase, Y (t), may take any value from the phase state space, Y, upon transitions [11]. The QBD process is called level-independent if all the transition rates (except the boundary) do not depend on X(t). The two-dimensional structure of the process allows to enumerate the states lexicographically into pairs (x, y), where x is the level value, and y is the phase value. Using this enumeration, the infinitesimal generator matrix, say Q, can be represented in the block-tridiagonal form:

Q =          A 0,0 A 0,1 0 0 . . . A 1,0 A (1) A (0) 0 . . . 0 A (2) A (1) A (0) . . . 0 0 A (2) A (1) . . . 0 0 . . . . . . . . .         

, where the element Q[(x 0 , y 0 ), (x 1 , y 1 )] of the matrix is the transition rate of the process {X(t), Y (t)} t 0 from the state (x 0 , y 0 ) to the state (x 1 , y 1 ) within adjacent levels,

|x 0 − x 1 | 1.

As such, A 0,0 is a square matrix consisting of transition rates at the (boundary) level 0, while A 0,1 (and A 1,0 , respectively) are the (possibly rectangular) matrices of outgoing (incoming) transitions from (to) level 0. The matrices A (2) , A (1) , A (0) are square matrices consisting of transition rates from level x to x − 1, x and x + 1, respectively. It is worth noting that only the diagonal elements of Q are negative, and

(A (0) + A (1) + A (2) )1 = 0.(1)

(In what follows, bold letter, say x, denotes a vector, x k is kth element of the vector, 1 is a vector of ones, and 0 is a vector of zeroes.) Commonly in the queueing applications, the level corresponds to the number of customers in the system or in the queue. Then transitions from level x to the level x − 1 and the level x + 1 correspond to a departure and an arrival, respectively. Being a special case of a continuous-time Markov chain, the irreducible QBD process has a specific form of the Foster negative drift criterion known as the Neuts ergodicity condition:

αA (0) 1 < αA (2) 1,(2)

where α 0 is the unique stochastic vector satisfying the system of linear algebraic equations α(A (0) + A (1) + A (2) ) = 0.

(

)3

Thus the QBD process is essentially considered at high levels when its stability is investigated.

In what follows, we consider the stability condition (2) of a level-independent irreducible QBD with the finite phase space Y. In general, the vector α in (3) has to be found numerically whenever the matrices A (0) , A (1) and A (2) are well defined. However, in the present paper the restrictions imposed on the possible phase transitions allow one to obtain α explicitly, even in the case when the matrix A (2) is not completely defined. Indeed, let A (0) and A (1) be diagonal (in such a case, the arrival process is a Poisson process with possibly phase-dependent rate). Thus at high levels the phase, Y (t), may change only at departure epochs.

Denote D = diag(d),

where d = A (2) 1 is the nonnegative vector of row sums of A (2) , i.e. d y is the total departure rate from phase y ∈ Y, and diag(d) is the square matrix consisting of zeroes except the main diagonal that contains the vector d. It then follows from (1) that A (1) = −A (0) − D.

Thus (3) may be written as αA (2) = αD,

which means that α depends only on the matrix A (2) . Intuitively this is expected, since the phase may change only at departure epochs. Moreover, it can be seen from ( 4) that the matrix A (2) does not necessarily need to be given explicitly. Indeed, l.h.s. of the first equation in ( 4) is a weighted column sum, whereas r.h.s. uses the row sums of the matrix A (2) . Define the matrix P = D −1 A (2) . Note that since A (0) and A (1) are diagonal, for the irreducible QBD the matrix D −1 exists. By noting that

D −1 y,y = d −1 y , y ∈ Y,

it can be easily seen that P is stochastic, since P 1 = 1. The matrix P contains probabilities of phase transitions of the Markov chain embedded at departure epochs. Denote by π the steady-state probability vector of P , that is, the solution of πP = π and π1 = 1. Let us introduce two new quantities:

α = CπD −1 , C = (πD −1 1) −1 =   y∈Y π y d −1 y   −1 .(5)

The expression for C follows from the last, balance equation of ( 4). It is easy to check by definition of D, P and π that α defined in ( 5) satisfies ( 4). The key observation that follows from ( 5) is that α depends only on the steady-state vector π and on the total departure rate vector d, and thus detailed knowledge of all elements of the matrix A (2) (or P ) is not necessary. Note that in such a case, the stability criterion (2) turns out to be

αA (0) 1 < C,(6)

where C is given by ( 5) and depends only on the vectors π and d. By rewriting (6) using α from (5), recalling that both D −1 and A (0) are diagonal, we obtain the following intuitively clear stability criterion:

y∈Y π y ρ y < 1,(7)

where

ρ y = A (0) y,y d y , y ∈ Y.

Indeed, ( 7) is a classical stability condition on the average load in the system that should be strictly less than one. In case the arrival process is Poisson, A (0) = λI (where I is identity matrix) and ( 6) transforms into

λ < C.(8)

However, the vector π has to be obtained given the matrix A (2) is not completely defined. In what follows, we obtain sufficient conditions for the model that allow π to be guessed.

Stability of a System with Automata-type Transitions

In this paper we consider the special type of queueing systems which we call the systems with automata-type transitions. Consider a queueing system in which customers of K classes arrive, 1 K < ∞. Customer's class is determined with a discrete distribution p = {p 1 , . . . , p K } at arrival epochs, which are time epochs of a Poisson process with possibly phase-dependent rate λ y , y ∈ Y. Classk customer requires exponentially distributed with rate µ k service time. The phase state space, Y, is the M -ary Cartesian power of the set {1, . . . , K}, i.e. Y ⊆ {1, . . . , K} M , and the phase y ∈ Y records the classes of M oldest (in the order of their arrival) customers in the system at time t ≥ 0 (from now on, in notation we stress that y is a vector). Customers never change their classes, while being among these M oldest. The classes of all customers being served at time t ≥ 0 are described in Y (t) in such a way that the service status of each of M oldest customers is unambiguously determined by the phase. More formally, there is an injection

σ : Y → {0, 1} M ,

such that σ m (y) = 1 if mth customer is being served if the system is in phase y, and σ m (y) = 0 if the mth customer is in the queue, for m = 1, . . . , M .

A particular example of such a system with automata-type transitions is the system with M servers in which each customer occupies at least one server. Here the information on classes of the M oldest customers present in the system (in the order of their arrival) and the total number in the system is sufficient to determine the next system state. Note that the classes of more recent customers (if any) are irrelevant for the system evolution. On the contrast, the two-class single-server system with absolute priorities is the example when the information on all the classes of customers present in the system is required, and thus such a system is not the one with the automata-type transitions and therefore is not in the focus of this paper.

The evolution of the considered system content can be described by the QBD process {X(t), Y (t)} t 0 , where X(t) is the total number of customers in the system, Y (t) = (Y 1 (t), . . . , Y M (t)) the phase that contains classes of the M oldest customers in the system at time t. In what follows, the M -dimensional vector y = (y 1 , . . . , y M ) is used to indicate the phase of the QBD in a lexicographical order, counting from the oldest customer. This QBD process has the following properties at high levels X(t) > M :

1. The phase is changed only at departure epochs. 2. At a departure epoch, the phase, y, evolves in such a way that only one component, say m, (corresponding to the departing customer) of y is removed (note that this means σ m (y) = 1), and the last component (corresponding to the class of M th oldest customer after transition) is being sampled from a discrete distribution p. Moreover, all the states in Y are communicating.

3. At a departure epoch, the total rate from state (x, y) is M m=1 µ ym σ m (y). These properties significantly restrict the structure of the generator submatrices A (i) , i = 0, 1, 2. Firstly, the matrices A (0) and A (1) must be diagonal (since the phase is changed only at departure epochs). Secondly, for y ∈ Y, the departure rate d y is given by

d y = M m=1 µ ym σ m (y).

Thirdly, the matrix A (2) is nonnegative with at least one positive element in each row.

Consider now a transition at a departure epoch, governed by the matrix P = D −1 A (2) . Fix y ∈ Y as the phase right after the transition. Recall that y M is the class of the M th oldest customer in the system after the transition, and this value is sampled from p. Denote by Y − (y) the set of states leading to y by a one-step transition, i.e. for any Y − (y) = {u ∈ Y : P (u, y) > 0}. Break the set Y − (y) into the subsets

Y − (y) = K k=1 Y − k (y),

using, as the index, the class, k, of the departing customer that caused the transition to y. The following lemma shows that these subsets are disjoint.

Lemma 1. For any y ∈ Y it holds that Y − i (y) ∩ Y − j (y) = ∅, i = j.

Proof. For any y ∈ Y define the characteristic vector

n(y) = (n 1 (y), . . . , n K (y)), n k (y) = M m=1 I(y m = k), k = 1, . . . , K,(9)

where I(•) is the indicator function. Note that n k (y) is the number of class-k customers among the M oldest in phase y.

From (9) it follows that

n(u) = n(k, y 1 , . . . , y M −1 ), u ∈ Y − k (y),

and

n k (u) 1. Denote n (k) = n(u), u ∈ Y − k (y). Then n(y) = n (k) − e k + e y M ,

where e k is an M -dimensional vector with 1 at the kth position and zero elsewhere. Thus for i = j, n (i) − e i = n (j) − e j , and hence n (i) = n (j) .

Due to the Lemma 1 from the balance equations π = πP we have

π y = K k=1 u∈Y − k (y) π u P (u, y), y ∈ Y.(10)

Take π y in the product form:

π y = M m=1 p ym , y = (y 1 , . . . , y M ) ∈ Y.(11)

Note that any element (i.e. phase) in the set Y − k (y) contains classes y 1 , . . . , y M −1 and the class of departing customer. Thus for any u ∈ Y − k (y), according to (11),

π u = p k M −1 m=1 p ym ,

and (10) transforms into

π y = K k=1 p k M −1 m=1 p ym   u∈Y − k (y) P (u, y)   .

Thus, if for all k = 1, . . . , K, y ∈ Y,

u∈Y − k (y) P (u, y) = p y M ,(12)

then the product form (11) is the solution of the balance equations.

In the next lemma we give the sufficient condition for (12) to hold.

µ k σ m (ϕ (m,k) (y)) d ϕ (m,k) (y) = 1. (13)

Proof. Fix y ∈ Y and k ∈ {1, . . . , K}. For any u ∈ Y − k (y) there exists such i ∈ {1, . . . , M } that u = ϕ (i,k) (y) and σ i (u) = 1 (possibly there are several such indices). Similarly, for a fixed i, by construction,

ϕ (i,k) (y) ∈ Y − k (y) if and only if σ i (ϕ (i,k) (y)) = 1. Thus, u∈Y − k (y) P u,y = M m=1 P (m) ϕ (m,k) (y),y ,(14)

where

P (m)

ϕ (m,k) (y),y is the probability of transition from the state ϕ (m,k) (y) to y by departure of mth oldest customer. The latter equals, by construction of the considered embedded Markov chain,

P (m) ϕ (m,k) (y),y = p y M µ k σ m (ϕ (m,k) (y)) d ϕ (m,k) (y) ,

and thus (13) together with (14) give (12). Now we formulate a sufficient condition for Lemma 2 to hold.

M m=1 σ m (ϕ (m,k) (y)) = σ k,y .(15)

Proof. Under the assumptions made,

d ϕ (m,k) (y) = µ k σ k,y ,

and hence

M m=1 µ k σ m (ϕ (m,k) (y)) d ϕ (m,k) (y) = M m=1 σ m (ϕ (m,k) (y)) σ k,y = 1.

4 Model Examples

Oldest Customer Class Chunk Policy

Consider the M -server system with K classes of customers, Poisson arrival process of rate λ, class-dependent service rates µ k and the following service policy: a group of the oldest customers of the same class are served, one customer per server. We may assume that a customer's class is sampled at M th position in the system (ordered by arrival times). In such a case, Y = {1, . . . , K} M is the set of classes of M oldest customers in the system. Fix some y ∈ Y and some class k. Then Y − k (y) contains only one state, u = (k, y 1 , . . . , y M −1 ). Moreover, σ(ϕ (m,k) (y)) = (1, 0, . . . , 0

) if k = y 1 , otherwise M m=1 σ m (ϕ (m,k) (y)) = M m=1 σ m (u) = max{i ≤ M : y 1 = • • • = y i = k} + 1.

Thus the stability criterion follows, after some algebra, from (7):

λ K k=1 1 µ k (1 − p k ) M −1 m=1 p m k m + p M k M < 1.

Weird Policy

Consider a two-server system with K = 2 types of customers. For y = (y 1 , y 2 ) ∈ Y = {1, 2} 2 define the following service rule:

σ(1, 1) = (0, 1), σ(1, 2) = (1, 0), σ(2, 1) = (1, 0), σ(2, 2) = (1, 1

). Thus, similar to the model in Section 4.1, for any y ∈ Y and any k = 1, 2, the set Y − k (y) contains only one state. Finally, the conditions of Corollary 1 are satisfied (which can be checked for each state y ∈ Y) and the product form stability criterion can be stated as

λ p 1 µ 1 + 2p 2 − p 2 2 2µ 2 < 1.

Supercomputer Model with Phase-Dependent Arrival Rate

In this section we reformulate the result first obtained in [17]. Consider an Mserver sytem with K = M classes of customers, where the input follows Poisson process of rate λ, and a customer of class k (arriving with probability p k ) requires k servers simultaneously for a random time, exponentially distributed with rate µ. The servers dedicated to a customer all are seized and released simultaneously. Clearly the service discipline is not work-conserving (a customer may be waiting in the queue when there are idle servers). As it was shown in [12,17], the system evolution may be described by the QBD process with X(t) being the number of customers in the system and Y (t) containing the classes of M oldest customers in the system (placed in the order of arrival). Thus for stability we consider the case X(t) > M and following the Corollary 1 check the sufficient condition (15) for the product form of π y . Fix y ∈ Y = {1, . . . , M } M and define σ m (y) = 1, m m(y), where

m(y) = max    i ≤ M : i j=1 y j ≤ M    . Fix k ≤ M and observe that for any u ∈ Y − k (y), M m=1 σ m (u) = σ k,y

and is constant. Indeed, for any such u,

m(u) = 1 + max    i ≤ M : k + i j=1 y j ≤ M    ,(16)

and hence m(u) depends only on k and y. It can be seen that m(u) = σ k,y . Finally, it follows from ( 16) that for any m ≤ M the value σ m (ϕ (m,k) (y)) = 1 if and only m ≤ σ k,y , and thus (15) follows.

Interestingly, the stability criterion of the supercomputer model in [17] was given in the form (8), where, as it was shown in [16], the value C may be obtained, after some algebra, in the form

C =   M m=1 1 mµ M j=m p * m j M k=M −j+1 p k   −1

, where the summation is performed over the number m of customers served, the number j of servers occupied by m customers, and the number k of servers required by the first customer waiting in the queue, and p * m j is the jth component of m times discrete convolution of the vector p with itself. However, the form of the result (7) allows to generalize the criterion to the model with phase-dependent input rate.

Indeed, let λ y be the input rate in the system if y ∈ Y is the phase of the system (we avoid the discussion of the service rate for the case X(t) < M since it is not used in the stability analysis, for completeness we may fix it as λ). Then A (0) y,y = λ y . It follows from (7) that the following stability criterion holds for the system with phase-dependent input rate:

y∈Y M m=1 p ym λ y m(y)µ < 1,

which is reduced to (8) if λ y ≡ λ.

Supercomputer Model with Different Service Rates

In this section we further generalize the supercomputer model by considering the case when classes have different service rates. It will be clear from subsequent analysis that the stability condition can be explicitly obtained for the model with M = 2 servers/classes. In this regard, we mention the pioneering work [6] where the 2-server model of this kind was first studied in a rather general context, for the general distribution of the class-2 customers (although the stability condition of such a case was not discussed extensively). Let customer of class m occupy m servers for exponentially distributed service time of rate µ m , m = 1, 2. Then the set Y contains only 4 phases, namely, the pairs (y 1 , y 2 ), y 1 , y 2 ∈ {1, 2}. It is also clear that σ(y) = (1, 1) only if y = (1, 1), and otherwise σ(y) = (1, 0). The subsequent analysis does not differ from Section 4.3, we only note that additionally we need to check the precondition of Corollary 1 that for any y ∈ Y and k ≤ 2, the rate µ um = µ k , for u ∈ Y − k (y) and m = 1, 2 such that σ m (u) = 1. This is straightforward, since the only phase where both customers are served is (1, 1), and the both have service rate µ 1 . However, this can not be generalized towards the system with M > 2 servers (if all classes are non-degenerate). It then follows that the stability condition of such a model with M = 2 servers is

p 2 1 λ (1,1) 2µ 1 + p 1 p 2 λ (1,2) µ 1 + p 1 p 2 λ (2,1) µ 2 + p 2 2 λ (2,2) µ 2 < 1,

which becomes λ < 2µ/(2 − p 2 1 ) (firstly appeared in [6]), if λ y ≡ λ and µ m ≡ µ, y ∈ Y, m = 1, 2.

Discussion

In the present paper we have introduced the special type of queueing systems with automata-type transitions. The idea is inspired by the concept of natural numbers and finite automata [2,3], where, loosely speaking, the output of a finite automation under certain conditions may have similar stochastic properties as the input. This theory was further developed into a concept of Stochastic Automata Networks [10] which also possess the product form solution under certain conditions on the automata interaction. We note that product form in many Markovian frameworks is tightly related to reversibility, either in time [5,18,21] or in structure [14]. In the present context, the reversibility may be thought of as time reversal, if for a state y we first sample the class of departing customer as k, and then sample the previous state as u ∈ Y − k (y) (note that the conditions of the Corollary 1 give the probabilities for such a sampling). This, however, requires proving that the probability of a departing (which corresponds to arrival backward in time) customer class k is p k , which we leave beyond the scope of this paper.

A natural question arises: is there an extension of the automata-type policy? An essential component is the possibility to make backward-in-time transitions, i.e. to define for any y ∈ Y the sets of states Y − k (y). This may be hard or even impossible, if the customers are ordered not in the order of arrival, or in any other deterministic way.

Another point of discussion is the connection between the QBD model and a more general simulation framework. Namely, if Y (t) is the phase process, and X(t) is the level process, then one can consider the so-called Generalized Semi-Markov Process related to the process {X(t), Y (t)} t≥0 in case the transitions are not Markovian. In such a case, the values σ i (y)µ yi may be thought of as clock speeds, and we need to keep track of the remaining work (which is not exponential in general), and possibly track the remaining interarrival time. Moreover, intuitively, the row-wise transformation of A (2) to P by multiplying on D −1 is inspired by the connection between time-and event-stationarity in a MAM framework, see e.g. [4]. At the same time, this is a more general property related to Palm theory, see [7,19]. This shows a direction for further study of the applicability of the obtained results to stochastic models beyond exponential distributions.

The possible generalizations of the model come from the point that at arrivals the phase is not changed. Thus, it seems technically straightforward to generalize the result towards MAP arrival processes (or even batch MAP arrivals), since the arrival phases will be unrelated to the phase of the system, and the rest will be done with Kronecker algebra. This is a promising direction of future research.