=Paper=
{{Paper
|id=Vol-2792/paper16
|storemode=property
|title=Stability of Multiclass Multiserver Models with Automata-type Phase Transitions
|pdfUrl=https://ceur-ws.org/Vol-2792/paper16.pdf
|volume=Vol-2792
|authors=Alexander Rumyantsev
}}
==Stability of Multiclass Multiserver Models with Automata-type Phase Transitions==
Stability of Multiclass Multiserver Models with
Automata-type Phase Transitions
Alexander Rumyantsev1,2
1
Institute of Applied Mathematical Research, Karelian Research Centre of RAS
2
Petrozavodsk State University
Tel.: +7-8142-76-63-12
ar0@krc.karelia.ru
Abstract. In this paper we obtain sufficient conditions for the multi-
class multiserver model with automata-type transitions to have a product-
form type explicit form of the stability criterion. We use the method to
obtain stability conditions of an interesting simultaneous service multi-
server model describing a supercomputer, both in case of phase-dependent
arrival rate and class-dependent service rate.
Keywords: Stochastic Stability · Multiclass Multiserver Model · Prod-
uct Form · Supercomputer Model · Simultaneous Service Multiserver
Model · Matrix-Analytic Method
1 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 two-
dimensional customers, namely, each such a customer requires a random number
of servers for a (same at each occupied server) random time, while servers ded-
icated to a customer are seized and released simultaneously. Such a model, to
Copyright © 2020 for this paper by its authors. Use permitted under
Creative Commons License Attribution 4.0 International (CC BY 4.0).
the best of our knowledge, for the first 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 flow (whereas the explicit form was given
for exponential service times, and for a particular case of phase-type distribu-
tion). In the present paper, as an example, we extend the results given in [17]
to a phase-dependent arrival flow and to a system with class-dependent service
times.
Summarizing, the contribution of this paper is two-fold. Firstly we obtain
such sufficient 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 sufficient 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 {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:
0,0 0,1
A A 0 0 ...
A1,0 A(1) A(0) 0 . . .
0 A(2) A(1) A(0) . . .
Q= ,
(1) . . .
0 (2)
0 A A
.. .. ..
0 0 . . .
where the element Q[(x0 , y0 ), (x1 , y1 )] of the matrix is the transition rate of
the process {X(t), Y (t)}t>0 from the state (x0 , y0 ) to the state (x1 , y1 ) within
adjacent levels, |x0 − x1 | 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)
(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 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 ma-
trix 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. dy 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, (4)
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
−1
Dy,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:
−1
X
α = CπD−1 , C = (πD−1 1)−1 = π y d−1
y
. (5)
y∈Y
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:
X
π y ρy < 1, (7)
y∈Y
where
(0)
Ay,y
ρy = , y ∈ Y.
dy
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 com-
pletely defined. In what follows, we obtain sufficient conditions for the model
that allow π to be guessed.
3 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 < ∞. Customer’s class is determined
with a discrete distribution p = {p1 , . . . , pK } at arrival epochs, which are time
epochs of a Poisson process with possibly phase-dependent rate λy , y ∈ Y. Class-
k 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) = (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 com-
ponent, 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.
PM
3. At a departure epoch, the total rate from state (x, y) is 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 dy is given by
M
X
dy = µym σ m (y).
m=1
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 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) = {u ∈ Y : P (u, y) > 0}. 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 ∈ Y it holds that Yi− (y) ∩ Yj− (y) = ∅, i 6= j.
Proof. For any y ∈ Y define the characteristic vector
M
X
n(y) = (n1 (y), . . . , nK (y)), nk (y) = I(ym = k), k = 1, . . . , K, (9)
m=1
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 ∈ Yk− (y),
and nk (u) > 1. Denote n(k) = n(u), u ∈ Yk− (y). Then
n(y) = n(k) − ek + eyM ,
where ek is an M -dimensional vector with 1 at the kth position and zero else-
where. Thus for i 6= j, n(i) − ei = n(j) − ej , and hence n(i) 6= n(j) . t
u
Due to the Lemma 1 from the balance equations π = πP we have
K
X X
πy = π u P (u, y), y ∈ Y. (10)
k=1 u∈Y − (y)
k
Take π y in the product form:
M
Y
πy = pym , y = (y1 , . . . , yM ) ∈ Y. (11)
m=1
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 ∈ Yk− (y), according to (11),
M
Y −1
π u = pk pym ,
m=1
and (10) transforms into
K
X M
Y −1 X
πy = pk pym P (u, y) .
k=1 m=1 u∈Yk− (y)
Thus, if for all k = 1, . . . , K, y ∈ Y,
X
P (u, y) = pyM , (12)
u∈Yk− (y)
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.
Lemma 2. For m = 1, . . . , M , k = 1, . . . , K and y ∈ Y construct the vectors
ϕ(m,k) (y) = (y1 , . . . , ym−1 , k, ym , . . . , yM −1 ). Then (12) holds if and only if for
any y ∈ Y
M
X µk σ m (ϕ(m,k) (y))
= 1. (13)
m=1
dϕ(m,k) (y)
Proof. Fix y ∈ Y and k ∈ {1, . . . , K}. For any u ∈ Yk− (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) ∈ Yk− (y) if and only
if σ i (ϕ(i,k) (y)) = 1. Thus,
M
(m)
X X
Pu,y = Pϕ(m,k) (y),y , (14)
u∈Yk− (y) m=1
(m)
where Pϕ(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,
(m) µk σ m (ϕ(m,k) (y))
Pϕ(m,k) (y),y = pyM ,
dϕ(m,k) (y)
and thus (13) together with (14) give (12). t
u
Now we formulate a sufficient condition for Lemma 2 to hold.
Corollary 1. Let for any y ∈ Y, k = 1, . . . , K, and u ∈ Yk− (y), µum = µk
PM
hold for any m = 1, . . . , M such that σ m (u) = 1. Let σ k,y := m=1 σ m (u) be
−
constant for all u ∈ Yk (y). Then (13) holds if
M
X
σ m (ϕ(m,k) (y)) = σ k,y . (15)
m=1
Proof. Under the assumptions made,
dϕ(m,k) (y) = µk σ k,y ,
and hence
M M
X µk σ m (ϕ(m,k) (y)) X σ m (ϕ(m,k) (y))
= = 1.
m=1
dϕ(m,k) (y) m=1
σ k,y
t
u
4 Model Examples
4.1 Oldest Customer Class Chunk Policy
Consider the M -server system with K classes of customers, Poisson arrival pro-
cess 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 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
X
σ m (ϕ(m,k) (y)) = σ m (u) = max{i ≤ M : y1 = · · · = yi = k} + 1.
m=1 m=1
Thus the stability criterion follows, after some algebra, from (7):
−1 m
K
" M
#
X 1 X pk pM
λ (1 − pk ) + k < 1.
µk m=1
m M
k=1
4.2 Weird Policy
Consider a two-server system with K = 2 types of customers. For y = (y1 , y2 ) ∈
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 Yk− (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
2p2 − p22
p1
λ + < 1.
µ1 2µ2
4.3 Supercomputer Model with Phase-Dependent Arrival Rate
In this section we reformulate the result first 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 sufficient condition (15)
for the product form of π y .
Fix y ∈ Y = {1, . . . , M }M and define σ m (y) = 1, m 6 m(y), where
X i
m(y) = max i ≤ M : yj ≤ M .
j=1
PM
Fix k ≤ M and observe that for any u ∈ Yk− (y), m=1 σ m (u) = σ k,y and is
constant. Indeed, for any such u,
Xi
m(u) = 1 + max i ≤ M : k + yj ≤ M , (16)
j=1
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
−1
M M M
X 1 X ∗m X
C= p pk ,
mµ j=m j
m=1 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 re-
quired 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
(0)
Ay,y = λy . It follows from (7) that the following stability criterion holds for the
system with phase-dependent input rate:
M
X Y λy
pym < 1,
m(y)µ
y∈Y m=1
which is reduced to (8) if λy ≡ λ.
4.4 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 (y1 , y2 ), y1 , y2 ∈ {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 ∈ 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
p21 λ(1,1) p1 p2 λ(1,2) p1 p2 λ(2,1) p2 λ(2,2)
+ + + 2 < 1,
2µ1 µ1 µ2 µ2
which becomes λ < 2µ/(2 − p21 ) (firstly appeared in [6]), if λy ≡ λ and µm ≡ µ,
y ∈ Y, m = 1, 2.
5 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 natu-
ral numbers and finite automata [2,3], where, loosely speaking, the output of a
finite automation under certain conditions may have similar stochastic proper-
ties as the input. This theory was further developed into a concept of Stochastic
Automata Networks [10] which also possess the product form solution under cer-
tain 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 ∈ Yk− (y) (note that the conditions
of the Corollary 1 give the probabilities for such a sampling). This, however, re-
quires 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 define for any y ∈ 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 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 exponen-
tial 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 re-
lated 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 19-
57-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.
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 finite automata. Soviet Mathematics Doklady
9, 324–325 (1968)
3. Becher, V., Heiber, P.A.: Normal numbers and finite 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. Brémaud, P.: Characteristics of queueing systems observed at events and the con-
nection 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 Au-
tomata 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 Can-
cer 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 Verification 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
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/0304-
4149(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 quasi-
reversible 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