Stability of Multiclass Multiserver Models with Automata-type Phase Transitions

Abstract

In this paper we obtain su cient conditions for the multiclass multiserver model with automata-type transitions to have a productform type explicit form of the stability criterion. We use the method to obtain stability conditions of an interesting simultaneous service multiserver model describing a supercomputer, both in case of phase-dependent arrival rate and class-dependent service rate.

Keywords

• Stochastic Stability uct Form
• Supercomputer Model
• Model
• Matrix-Analytic Method

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 speci c 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 the best of our knowledge, for the rst time was treated by MAM in [12], while in [6] the stability condition of a two-server system was expressed in explicit form, whereas in [8] the condition was expressed for a two-server system with Markovian arrival process (MAP). In [17] the stability criterion in explicit form was obtained for arbitrary number of servers, Poisson arrivals and exponentially distributed class-independent service times, while in [1] it was formulated in a general case for a regenerative input ow (whereas the explicit form was given for exponential service times, and for a particular case of phase-type distribution). In the present paper, as an example, we extend the results given in [17] to a phase-dependent arrival ow and to a system with class-dependent service times.

Summarizing, the contribution of this paper is two-fold. Firstly we obtain such su cient conditions for the multiclass multiserver model that the stability criterion can be obtained explicitly. Second, we extend the stability criterion of a simultaneous service multiserver model obtained in [17] towards phase-dependent arrival rate system and/or system with class-dependent service rate.

The structure of the paper is as follows. In Section 2 we review the stability criterion of a MAM known as the Neuts ergodicity condition. In Section 3 we give the description of a multiclass multiserver model with automata-type transitions and formulate the conditions for the stability condition to be given in product form. In Section 4 we demonstrate the application of the su cient condition obtained in Section 3 to various models including supercomputer model with phase-dependent arrival rate and/or class-dependent service rate, and discuss the limitations of such conditions. The discussion of potential extensions of the work is given in Section 5. 2

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 fX(t); Y (t)gt>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 in nitesimal generator matrix, say Q, can be represented in the block-tridiagonal form:

A( 2 ) A( 1 ) A(0) . . . CCC ;

C 0 A( 2 ) A( 1 ) . . . CC 0 . . . . . . . . . A 0 : : : 1 0 : : : C (In what follows, bold letter, say x, denotes a vector, xk 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 speci c form of the Foster negative drift criterion known as the Neuts ergodicity condition: where > 0 is the unique stochastic vector satisfying the system of linear algebraic equations

A(0)1 <

A( 2 )1; (A(0) + A( 1 ) + A( 2 )) = 0: where the element Q[(x0; y0); (x1; y1)] of the matrix is the transition rate of the process fX(t); Y (t)gt>0 from the state (x0; y0) to the state (x1; y1) within adjacent levels, jx0 x1j 6 1. As such, A0;0 is a square matrix consisting of transition rates at the (boundary) level 0, while A0;1 (and A1;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 ) ( 2 ) ( 3 ) ( 4 ) 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 nite 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 de ned. 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 de ned. 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. dy is the total departure rate from phase y 2 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

Thus (3) may be written as

A( 1 ) =

A(0)

D: A( 2 ) = 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 rst equation in ( 4 ) is a weighted column sum, whereas r.h.s. uses the row sums of the matrix A( 2 ).

De ne the matrix P = D 1A( 2 ). Note that since A(0) and A( 1 ) are diagonal, for the irreducible QBD the matrix D 1 exists. By noting that

Dy;1y = d 1 y ; y 2 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 11) 1 = 4X ydy 15

: 2 y2Y The expression for C follows from the last, balance equation of ( 4 ). It is easy to check by de nition of D, P and that de ned in ( 5 ) satis es ( 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 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: where

A(0)1 < C; X y2Y

y y < 1; y =

A(y0;y) ; dy

y 2 Y: < C: 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

However, the vector has to be obtained given the matrix A( 2 ) is not completely de ned. In what follows, we obtain su cient conditions for the model that allow to be guessed. ( 5 ) ( 6 ) ( 7 ) ( 8 )

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 6 K < 1. Customer's class is determined with a discrete distribution p = fp1; : : : ; pK g at arrival epochs, which are time epochs of a Poisson process with possibly phase-dependent rate y; y 2 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 f1; : : : ; Kg, i.e.

Y

f1; : : : ; KgM ; and the phase y 2 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 ! f0; 1gM ; 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 su cient 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 fX(t); Y (t)gt>0, where X(t) is the total number of customers in the system, Y (t) = (Y1(t); : : : ; YM (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 = (y1; : : : ; yM ) 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 PM m=1 ym m(y).

These properties signi cantly 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 2 Y, the departure rate dy is given by

M dy = X 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 1A( 2 ). Fix y 2 Y as the phase right after the transition. Recall that yM 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) = fu 2 Y : P (u; y) > 0g. Break the set Y (y) into the subsets

K Y (y) = [ Yk (y);

k=1 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 2 Y it holds that Yi (y) \ Yj (y) = ?; i 6= j. Proof. For any y 2 Y de ne the characteristic vector n(y) = (n1(y); : : : ; nK (y)); nk(y) = M X I(ym = k); k = 1; : : : ; K; m=1 ( 9 ) where I( ) is the indicator function. Note that nk(y) is the number of class-k customers among the M oldest in phase y.

From ( 9 ) it follows that

n(u) = n(k; y1; : : : ; yM 1); u 2 Yk (y); and nk(u) > 1. Denote n(k) = n(u); u 2 Yk (y). Then n(y) = n(k) ek + eyM ; where ek is an M -dimensional vector with 1 at the kth position and zero elsewhere. Thus for i 6= j, n(i) ei = n(j) ej, and hence n(i) 6= n(j). tu Due to the Lemma 1 from the balance equations = P we have

K y = X

X k=1 u2Yk (y) uP (u; y); y 2 Y: ( 10 )

Take

y in the product form:

M y = Y pym ; m=1 y = (y1; : : : ; yM ) 2 Y: ( 11 ) Note that any element (i.e. phase) in the set Yk (y) contains classes y1; : : : ; yM 1 and the class of departing customer. Thus for any u 2 Yk (y), according to ( 11 ), and ( 10 ) transforms into Thus, if for all k = 1; : : : ; K; y 2 Y; y =

K M 1 X pk Y k=1 m=1

3

P (u; y)5 :

M 1 u = pk Y pym ; m=1

2 pym 4

X u2Yk (y)

X u2Yk (y)

P (u; y) = pyM ; then the product form ( 11 ) is the solution of the balance equations.

In the next lemma we give the su cient condition for ( 12 ) to hold. Lemma 2. For m = 1; : : : ; M , k = 1; : : : ; K and y 2 Y construct the vectors '(m;k)(y) = (y1; : : : ; ym 1; k; ym; : : : ; yM 1). Then ( 12 ) holds if and only if for any y 2 Y

M X m=1 k m('(m;k)(y)) d'(m;k)(y) = 1: Proof. Fix y 2 Y and k 2 f1; : : : ; Kg. For any u 2 Yk (y) there exists such i 2 f1; : : : ; M g that u = '(i;k)(y) and i(u) = 1 (possibly there are several such indices). Similarly, for a xed i, by construction, '(i;k)(y) 2 Yk (y) if and only if i('(i;k)(y)) = 1. Thus,

X u2Yk (y)

M Pu;y = X m=1

P('m(m) ;k)(y);y; 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 = pyM and thus ( 13 ) together with ( 14 ) give ( 12 ).

Now we formulate a su cient condition for Lemma 2 to hold. ( 12 ) ( 13 ) ( 14 ) Corollary 1. Let for any y 2 Y, k = 1; : : : ; K, and u 2 Yk (y), hold for any m = 1; : : : ; M such that m(u) = 1. Let constant for all u 2 Yk (y). Then ( 13 ) holds if k;y := PM m=1

Proof. Under the assumptions made, and hence

4 4.1 M X m=1 m('(m;k)(y)) = k;y: d'(m;k)(y) = k k;y; M X m=1 ( 15 ) tu

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 = f1; : : : ; KgM is the set of classes of M oldest customers in the system. Fix some y 2 Y and some class k. Then Yk (y) contains only one state, u = (k; y1; : : : ; yM 1). Moreover, ('(m;k)(y)) = (1; 0; : : : ; 0) if k 6= y1, otherwise

M X m=1

M m('(m;k)(y)) = X m=1 m(u) = maxfi

M : y1 = = yi = kg + 1: Thus the stability criterion follows, after some algebra, from ( 7 ): XK 1 " k=1 k (1

M 1 pm pk) X

k + m=1 m pM # k M < 1: 4.2

Weird Policy

Consider a two-server system with K = 2 types of customers. For y = (y1; y2) 2 Y = f1; 2g2 de ne 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 2 Y and any k = 1; 2, the set Yk (y) contains only one state. Finally, the conditions of Corollary 1 are satis ed (which can be checked for each state y 2 Y) and the product form stability criterion can be stated as < 1: 4.3

Supercomputer Model with Phase-Dependent Arrival Rate In this section we reformulate the result rst obtained in [17]. Consider an M server sytem with K = M classes of customers, where the input follows Poisson process of rate , and a customer of class k (arriving with probability pk) 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 su cient condition ( 15 ) for the product form of y.

Fix y 2 Y = f1; : : : ; M gM and de ne m(y) = 1; m 6 m(y), where Fix k M and observe that for any u 2 Yk (y), PmM=1 constant. Indeed, for any such u, m(u) =

k;y and is 8 < m(u) = 1 + max i :

M : k +

i X yj j=1

M

; 9 = ; ( 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 = 42 XM m1 XM pjm m=1 j=m

M

X k=M j+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 rst customer waiting in the queue, and pjm 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 2 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 x it as ). Then A(y0;)y = y. It follows from ( 7 ) that the following stability criterion holds for the system with phase-dependent input rate:

M X Y In this section we further generalize the supercomputer model by considering the case when classes have di erent 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 rst 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 (y1; y2); y1; y2 2 f1; 2g. It is also clear that (y) = (1; 1) only if y = (1; 1), and otherwise (y) = (1; 0). The subsequent analysis does not di er from Section 4.3, we only note that additionally we need to check the precondition of Corollary 1 that for any y 2 Y and k 2, the rate um = k, for u 2 Yk (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 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 nite automata [2,3], where, loosely speaking, the output of a nite 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] which becomes y 2 Y; m = 1; 2. 5

Discussion < 2 =(2 p21) ( rstly appeared in [6]), if y and m or in structure [14]. In the present context, the reversibility may be thought of as time reversal, if for a state y we rst sample the class of departing customer as k, and then sample the previous state as u 2 Yk (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 pk, 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 de ne for any y 2 Y the sets of states Yk (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 SemiMarkov Process related to the process fX(t); Y (t)gt 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. Acknowledgements This work is partially supported by RFBR, projects 1957-45022, 19-07-00303. Author would like to thank Prof. Evsey Morozov, Prof. Miklos Telek, Prof. Srinivas R. Chakravarthy and anonymous referees for helpful discussions and suggestions that improved the paper. 18. Taylor, P.: Quasi-reversibility and networks of queues with nonstandard batch movements. Mathematical and Computer Modelling 31( 10-12 ), 335{ 341 (May 2000). https://doi.org/10.1016/S0895-7177(00)00104-7, https:// linkinghub.elsevier.com/retrieve/pii/S0895717700001047 19. Thorisson, H.: On time- and cycle-stationarity. Stochastic Processes and their Applications 55( 2 ), 183{209 (Feb 1995). https://doi.org/10.1016/03044149(94)00038-U, http://linkinghub.elsevier.com/retrieve/pii/ 030441499400038U 20. Tweedie, R.L.: Relations between ergodicity and mean drift for markov chains. Australian Journal of Statistics 17( 2 ), 96{102 (Jun 1975). https://doi.org/10.1111/j.1467-842X.1975.tb00945.x, https://doi.org/10. 1111/j.1467-842X.1975.tb00945.x, publisher: John Wiley & Sons, Ltd 21. Walrand, J., Varaiya, P.: Interconnections of Markov chains and quasireversible queuing networks. Stochastic Processes and their Applications 10( 2 ), 209{219 (Sep 1980). https://doi.org/10.1016/0304-4149(80)90022-8, http://www. sciencedirect.com/science/article/pii/0304414980900228

References

1. Afanaseva , L. , Bashtova , E. , Grishunina , S. : Stability Analysis of a Multi-server Model with Simultaneous Service and a Regenerative Input Flow . Methodology and Computing in Applied Probability (Jul 2019 ). https://doi.org/10/ggnsgj, http: //link.springer. com/10.1007/s11009-019-09721-9

2. Agafonov , V.: Normal sequences and nite automata . Soviet Mathematics Doklady 9 , 324 { 325 ( 1968 )

3. Becher , V. , Heiber , P.A. : Normal numbers and nite automata . Theoretical Computer Science 477 , 109 {116 (Mar 2013 ). https://doi.org/10/f4rn5f, https: //linkinghub.elsevier.com/retrieve/pii/S0304397513000698

4. Bladt , M. , Nielsen , B.F. : Matrix-Exponential Distributions in Applied Probability, Probability Theory and Stochastic Modelling , vol. 81 . Springer US , Boston, MA ( 2017 ). https://doi.org/10.1007/978-1- 4939 -7049-0, http://link.springer.com/ 10.1007/978-1- 4939 -7049-0

5. Brill , P.H. , Cheung , C. h ., Hlynka , M. , Jiang , Q. : Reversibility Checking for Markov Chains . Communications on Stochastic Analysis 12 ( 2 ) (Oct 2018 ). https://doi.org/10.31390/cosa.12.2.02, https://digitalcommons.lsu.edu/cosa/ vol12/iss2/2

6. Brill , P.H. , Green , L. : Queues in which customers receive simultaneous service from a random number of servers: a system point approach . Management Science 30 ( 1 ), 51 { 68 ( 1984 ). https://doi.org/10.1287/mnsc.30.1.51, http://pubsonline. informs.org/doi/abs/10.1287/mnsc.30.1. 51

7. Bremaud , P. : Characteristics of queueing systems observed at events and the connection between stochastic intensity and Palm probability . Queueing systems 5(1- 3) , 99 { 111 ( 1989 ). https://doi.org/10.1007/BF01149188, http://link.springer. com/article/10.1007/BF01149188

8. Chakravarthy , S.R. , Karatza , H.D.: Two-server parallel system with pure space sharing and Markovian arrivals . Computers & Operations Research 40 ( 1 ), 510 { 519 ( 2013 ). https://doi.org/10.1016/j.cor. 2012 . 08 .002

9. Foster , F.G. : On the Stochastic Matrices Associated with Certain Queuing Processes . Ann. Math. Statist. 24 ( 3 ), 355 {360 (Sep 1953 ). https://doi.org/10.1214/aoms/1177728976, https://projecteuclid.org: 443/euclid.aoms/1177728976, publisher: The Institute of Mathematical Statistics

10. Fourneau , J. , Plateau , B. , Stewart , W. : Product form for Stochastic Automata Networks . In: Proceedings of the 2nd International ICST Conference on Performance Evaluation Methodologies and Tools . ICST, Nantes, France ( 2007 ). https://doi.org/10.4108/valuetools. 2007 . 1980 , http://eudl.eu/doi/10. 4108/valuetools. 2007 .1980

11. He , Q.M.: Fundamentals of Matrix-Analytic Methods . Springer New York ( 2014 )

12. Kim , S.: M/M/s Queueing System Where Customers Demand Multiple Server Use . Ph.D. thesis , Southern Methodist University ( 1979 )

13. Meyn , S.P. , Tweedie , R.L. : Markov chains and stochastic stability . Springer Science & Business Media ( 2012 )

14. Miyazawa , M.: Reversibility in Queueing Models . In: Wiley Encyclopedia of Operations Research and Management Science , pp. 1 { 19 . American Cancer Society ( 2013 ). https://doi.org/10.1002/9780470400531.eorms1079, https:// onlinelibrary.wiley.com/doi/abs/10.1002/9780470400531.eorms1079

15. Neuts , M.F. : Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Approach . Johns Hopkins University Press, Baltimore ( 1981 )

16. Rumyantsev , A. , Morozov , E.: Accelerated Veri cation of Stability of Simultaneous Service Multiserver Systems . In: Proceedings of 2015 7th International Congress on Ultra Modern Telecommunications and Control Systems and Workshops (ICUMT) , Brno, Czech Republic, 6 - 8 October 2015 . pp. 239 { 242 . IEEE ( 2015 )

17. Rumyantsev , A. , Morozov , E.: Stability criterion of a multiserver model with simultaneous service . Annals of Operations Research 252 ( 1 ), 29 { 39 ( 2017 ). https://doi.org/10.1007/s10479-015-1917-2, https://doi.org/10.1007/ s10479-015-1917-2