Systems Engineering - Theory & Practice Volume 27, Issue 1, January 2007 Online English edition of the Chinese language journal Cite this article as: SETP, 2007, 27(1): 91–98
− + Geo2/Geo1, Geo2/s/s + K Queue System with Mixed Geo1← Delay and Loss MA Zhan-you1 , 2 ,∗ , LIU Ming-xin3 , XU Xiu-li1 , TIAN Nai-shuo1 1. College of Sciences, Yanshan University, Qinhuangdao 066004, China 2. College of Mathematics, Jilin Normal University, Siping 136000, China 3. College of Information Science and Engineering, Yanshan University, Qinhuangdao 066004, China −Geo /Geo , Geo /s/s + K queue system with two classes of customers and priorities is considered. Abstract: In this study, a Geo1 ← + 2 1 2
The first class of customers are with delay and loss, and the second class of customers are with loss. The second class of customers have preemptive resume priority to the first class of customers. Using the method of matrix analysis, the mean queue length, the loss probability and the system utility for the two classes of customers are derived. Finally, some numerical results are presented by Matlab program. Key Words: discrete time queue; preemptive resume priority; mixed policy; immediate access; matrix-geometric solution
1 Introduction Unitary standard of QoS cannot satisfy the requirements for some customers because of the diversity of quality of service (QoS) for different kinds of customers; therefore, it impels us to classify customers by priorities and to provide different QoS for different classes of customers based on their priorities. The priority queues are considered to be more appropriate for modeling practical problems than others, such as, communication networks and sales of perishables. Douglas[1] studied an M/M/1 queue with single server and infinity waiting space to describe a priority model with two classes of customers. Andreas et al[2] and Hu[3] discussed an M/M/s priority queue with two classes of customers, Niu[4] investigated a continuous time-priority queue with two classes of customers and infinity waiting space, and also got stationary performance indices of the system. However, majority of literatures only analyzed continuous time priority queues, for example, literatures [5–6]; moreover, cell arrivals take place in discrete instants with fixed interarrivals in communication networks. Therefore, this study examines a discrete time multiserver queue with two classes of customers and preemptive resume priority using matrixgeometric solution method[7−10] .
P {S = k} = µ ¯k−1 µi , i = 1, 2, k = 1, 2, · · ·. The arrivals in i different slots and the arrivals of different class of customers in the same slot are independent. The priority policy of the system is shown in Figure 1. The second class of customers have the highest priority, that is, if a second class of customers arrives in the system and all servers are busy and parts of servers are serving the first class of customers, one of the servers who serves the first class of customers completes serving immediately and serves this newly arrived customer; if a second class of customers arrives in the system and all servers are busy serving the second class of customers, this newly arrived customer is rejected. When all servers are busy serving the first class of customers, a customer of the first class who arrives or who is interrupted serving by the second class of customers has to queue in the limit waiting space. The model is referred to as preemptive priority queue with loss and delay (PP-L-D), −Geo /Geo , Geo /s/s + K. + denoted by Geo1 ← 2 1 2
2 Model description Assume that there are s servers to provide service for two classes of customers in the system. The first class of customers are with delay and loss (for example, date) and the second class of customers are with loss (for example, sound, image). The interarrival times and service times of the ith class of customers follow geometric distributions with pa¯ k−1 λi , rameter λi and µi , respectively, i.e. P {T = k} = λ i
Figure 1. Priority policy of the system
On the basis of the states of the second class of cus-
Received date: July 8, 2005 ∗ Corresponding author: Tel: +86-335-8074703; E-mail:
[email protected] Foundation item: Supported by the National Natural Science Foundation of China (No.10271102, 10671170) c 2007, Systems Engineering Society of China. Published by Elsevier BV. All rights reserved. Copyright
MA Zhan-you, et al./Systems Engineering – Theory & Practice, 2007, 27(1): 91–98
tomers, there are three cases in the system: State I (idle state): Parts or all s servers are idle and a new arrival can be served immediately; State II (priority state): s servers are all busy and the number of the second class of customers is less than s in the system. A new arrival of the second class preempts the server, which is serving for the first class of customers; hence, this customer whose service is interrupted has no choice but to wait in the front of the queue. A new arrival of the first class also has to wait in line. State III (barrage state): s servers are all serving the second class of customers, a new arrival of the second class is barraged, and a new arrival of the first class has to queue in the waiting space.
3 Transition probability matrix P The time axis is divided into small section with fixed length, which is called a slot. The arrivals of customers only take place at the end of every slot t = n− , n = 1, 2, · · ·, this model is referred to as the late-arrival system. To define the number of customers in the discrete instant t = n in the system distinctly, we assume that the beginning and ending of service happen in the instant t = n. The customer who arrived in the time n− can be served in the time t = n immediately if a customer has departed in the time t = n. This model is referred to as a system with immediate access (see Hunter[10] ). Assume that interarrival times and service times are independent, let Ln = L(n+ ) be the number of customers (the sum of two classes of customers) in the instant t = n+ in the (i) system, and Ln = L(i) (n+ ) be the number of the ith class of customers in the instant t = n+ in the system, i = 1, 2. According to arrivals and departures of the customers, the customer arrived in the time t = n− can be reckoned in Ln , and the customer departed in the time t = n does not so. If (2) {(Ln , Ln ), n ≥ 1} is a bidimensional Markov chain with state space
then the transition probability matrix of the system can be written as P0,0 P0,1 P0,2 ⎜ .. .. .. .. ⎜. . . . ⎜ ⎜ Ps,0 Ps,1 · · · Ps,s+2 P =⎜ ⎜ Ps+1,1 Ps+1,2 · · · ⎜ ⎜ .. .. ⎝ . .
⎞
Ps+1,s+3 .. .
Ps+K,K−1 · · ·
Pi,i+2
¯i1 λ2 µ ¯02 0 λ1 µ ⎜ .. = ⎝ ... . 0 λ1 µ ¯01 λ2 µ ¯i2
⎞ 0 .. ⎟ . ⎠ 0
(i+1)×(i+3)
if i = s − 1, we have ⎛ ¯s−1 λ2 µ ¯02 0 λ1 µ 1 ⎜ .. .. Ps−1,s+1 = ⎝ . .
⎞ ⎟ ⎠ ¯01 λ2 µ ¯s−1 0 λ1 µ 2
s×(s+1)
if i = s, we obtain ⎞ ⎛ ¯s1 λ2 µ ¯02 0 λ1 µ ⎟ ⎜ .. .. ⎟ ⎜ . . Ps,s+2 = ⎜ ⎟ ⎠ ⎝ ¯ 1 λ2 µ ¯s−1 0 λ1 µ 2 0 (s+1)×(s+1) (2) If 0 ≤ i ≤ s − 1, Pi,i+1 are (i + 1) × (i + 2) matrix, furthermore, Ps,s+1 is (s + 1) × (s + 1) matrix, we obtain ⎞ ⎛ 0 1 1 0 Bi D0 Bi D0 ⎟ ⎜ .. .. Pi,i+1 = ⎝ ⎠ . . B00 Di1 B01 Di0
(i+1)×(i+2)
⎞ Bs0 D01 Bs1 D00 ⎟ ⎜ .. .. ⎟ ⎜ . . =⎜ ⎟ 0 1 1 0 ⎠ ⎝ B1 Ds−1 B1 Ds−1 0 1 0 B0 (Ds + Ds ) (s+1)×(s+1) ⎛
Ps,s+1
(3) If 0 ≤ i ≤ s, Pi,0 is (i + 1) × 1 matrix, then their transposed matrix (Pi,0 )T are as follows ¯ 1 µi λ ¯ 0¯ i ¯ 0 ¯ i−1 ¯ (Pi,0 )T= (λ 1 2 µ2 , λ1 µ1 λ2 µ2 , · · · , λ1 µ1 λ2 µ2 )1×(i+1)
Ω = {(n, k) : 0 ≤ n ≤ s − 1, 0 ≤ k ≤ n} ∪{(n, k) : s ≤ n ≤ s + K, 0 ≤ k ≤ s}.
⎛
⎛
⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ Ps+K,s+K (1)
j ¯i−j To express every subblock in P , denoted bji = Cij µ 1 µ1 , j j i−j j j j 0 0 ¯2 µ2 , especially, bi = di = 0, b0 = d0 = di = Ci µ 1, when j > i or j < 0, i ≥ 0. hence, we obtain ¯ 1 bj−1 = B j , λ2 dj + λ ¯ 2 dj−1 = Dj , where eleλ1 bji + λ i i i i i ments in P are denoted as follows: (1) If 0 ≤ i ≤ s − 2, Pi,i+2 are (i + 1) × (i + 3) matrix and can be written as
(4) If 2 ≤ i ≤ s − 1, Pi,i are (i + 1)-order tridiagonal matrix ,we have ¯ ¯ P0,0 = λ 1 λ2 ¯1µ ¯2 λ ¯ 1 µ1 λ2 ¯ 1 )λ (λ1 µ1 + λ P1,1 = ¯2µ ¯ ¯ 1 (λ2 µ2 + λ ¯2⎞) λ ⎛ 1 λ11λ2 µ2 2 0 Bi D0 Bi D 0 0 1 2 ⎟ ⎜Bi−1 D12 Bi−1 D11 Bi−1 D10 ⎟ ⎜ Pi,i=⎜ ⎟ .. .. .. ⎠ ⎝ . . . 0 2 1 1 B0 Di B0 Di (i+1)×(i+1) ⎛ 1 1 Bs D0 Bs2 D00 0 1 2 ⎜Bs−1 D12 Bs−1 D11 Bs−1 D10 ⎜ Ps,s=⎜ .. .. .. ⎝ . . . λ1 Ds2
⎞
⎟ ⎟ ⎟ ⎠ 0 1 0 B0 (Ds + Ds ) (s+1)×(s+1)
(5) If 2 ≤ i ≤ s, 1 ≤ j < i, Pi,j are (i + 1) × (j + 1) matrix, we have
MA Zhan-you, et al./Systems Engineering – Theory & Practice, 2007, 27(1): 91–98
⎛
Pi,j
Bii+1−j D01 ⎜ . ⎜ .. ⎜ ⎜ B i−2j+2 Dj ⎜ i−j+1 j−1 ⎜ . =⎜ ⎜ .. ⎜ 0 i−j+2 ⎜ Bj−1 Di−j+1 ⎜ ⎜ ⎝
Bii+2−j D00 .. .
i−2j+3 j−1 Bi−j+1 Dj−1 .. . i−j+1 1 Bj−1 Di−j+1
⎞ .. . i−2j+4 j−2 Bi−j+1 Dj−1 .. . i−j 2 Bj−1 Di−j+1 .. .
(6) If s + 1 ≤ i ≤ 2s, Pi,j is denoted by
Pi,j =
i−s
i−j+2 0 Bi−j+1 Dj−1 .. . j+1 i−2j+2 Bj−1 Di−j+1 .. . B01 Dii−j+1
⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (i+1)×(j+1)
Stationary probability distribution
From equilibrium equation ΠP = Π and normalization condition Πe = 1, we can easily obtain the stationary probability distribution under the equilibrium condition, where e is a column vector of ones. According to the relationship between K and s, we discuss the stationary probability distribution in three cases:
s
where (Ps,j+s−i , 0, · · · , 0), is composed of matrix Ps,j+s−i
i−s
and i − s column vectors of zeros. (7) If 2s < i < s + K − 1, we obtain P2s,j+2s−i , i − s ≤ j ≤ i + 2 Pi,j = 0, other
Theorem 4.1 If 1 ≤ K ≤ s, Π satisfies equilibrium equation ΠP = Π and normalization condition Πe = 1, then we have
(8) The subblock in the last two rows of the transition probability matrix can be written as ⎧ K −1≤j ≤s+K −1 ⎨ Ps+K−1,j , Ps,s+1 + Ps,s+2 , j = s + K Ps+K−1,j = ⎩ 0, other
Ps+K,j
. ··· .. . ··· .. . B00 Dii−j+2
4
⎧ (Ps,s−1 , 0), j =i−1 ⎪ ⎪ ⎪ ⎪ , 0, · · · , 0 ), i−s≤j ≤s (P ⎨ s,j+s−i
⎪ ⎪ ⎪ , P ⎪ ⎩ s,j+s−i 0,
..
⎧ K ≤j ≤s+K −1 ⎨ Ps+K,j , Ps,s + Ps,s+1 + Ps,s+2 , j = s + K = ⎩ 0, other
(π0 , π1 , · · · , πs+K )P = (π0 , π1 , · · · , πs+K ) (π0 , π1 , · · · , πs+K )e = 1
(2)
where e is a column vector of ones. The proof of Theorem 4.1 is similar to literature [6]. Lemma 4.1 If K > s, Markov chain {(L, L(2) ), n ≥ 1} of transition probability matrix P is positive recurrent if and only if SP (R) < 1, where SP (R) is spectrum radius of R, which is the smallest nonnegative solution of matrix equation s+2
Ri P2s,2s+2−i = 0.
(3)
i=0
The structure of the transition probability matrix P indicates that {(L(n), L(2) (n)), n ≥ 1} is nonperiodic and irreducible; therefore, {(L(n), L(2) (n)), n ≥ 1} is also positive recurrent. {(L, L(2) )} is the stationary limiting of bidimensional Markov chain {(L(n), L(2) (n)), n ≥ 1}, denoted Π = {π0 , π1 , · · · , πs , πs+1 , · · · , πs+K }, where every component is as follows πk = (πk,0 , πk,1 , · · · , πk,k ), πk = (πk,0 , πk,1 , · · · , πk,s ), ⎛
P0,0 ⎜ .. ⎜ . ⎜ ⎜ Ps,0 ⎜ B[R] = ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
Lemma 4.2 If matrix G = P2s,s + P2s,s+1 + · · · + P2s,2s+2 is irreducible, rate matrix R satisfies the condition SP (R) < 1 if and only if π ∗ (P2s,2s+1 + 2P2s,2s+2 )e < π ∗ (sP2s,s + (s − 1)P2s,s+1 + · · · + P s2s,2s−1 )e, where π ∗ is the stationary probability vector of stochastic matrix G, and satisfies equations π ∗ G = π ∗ , π ∗ e = 1. The proofs of Lemma 4.1 and Lemma 4.2 are given in literature [9].
0 ≤ k < s,
Theorem 4.2 If SP (R) < 1 and s < K ≤ 2s, stochastic matrix
s ≤ k ≤ s + K.
s(3s+1) -order 2
⎞
P0,1 .. .
P0,2 .. .
Ps,1 Ps+1,1
Ps,2 Ps+1,2 .. .
. ··· ··· .. .
P2s,s
···
..
Ps,s+1 Ps+1,s+1 K i=0
RK−i P2s,s+i+1
Ps,s+2 Ps+1,s+2 .. . ···
Ps+1,s+3 .. . K i=0
RK−i P2s,2s+i−K
⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
MA Zhan-you, et al./Systems Engineering – Theory & Practice, 2007, 27(1): 91–98
has positive vector (π0 , π1 , · · · , π2s ) such that (π0 , π1 , · · · , π2s )B[R] = (π0 , π1 , · · · , π2s )
(4)
simultaneously, the stationary probability distribution can be denoted as πk = π2s Rk−2s ,
2s < k ≤ s + K − 1
and (π0 , π1 , · · · , π2s ) is determined by the following set of equations ⎛
P0,0 ⎜ .. ⎜ . ⎜ ⎜ Ps,0 ⎜ B[R] = ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
P0,1 .. .
P0,2 .. .
Ps,1 Ps+1,1
Ps,2 Ps+1,2 .. .
⎧ (π0 , π1 , · · · , π2s )B[R] = (π0 , π1 , · · · , π2s ) ⎪ ⎪ ⎪ 2s−1 ⎪ ⎨ πi e + π2s (I − RK−s )(I − R)−1 + πs+K e = 1 i=0 ⎪ ⎪ = (πs+K−1 Ps+K−1,s+K + πs+K−2 Ps+K−2,s+K ) π ⎪ ⎪ ⎩ s+K (I − Ps+K,s+K )−1 (5) s(3s+1) 2
Theorem 4.3 If SP (R) < 1 and K > 2s, order stochastic matrix ⎞
··· ··· .. .
Ps,s+2 Ps+1,s+2
P2s,s
···
Ps+1,s+3 .. . s−1 i=0
Rs−1−i P2s,2s−1
s i=0
Rs−i P2s,s+i
⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
has positive vector (π0 , π1 , · · · , π2s ) such that (π0 , π1 , · · · , π2s )B[R] = (π0 , π1 , · · · , π2s )
P (L(1) =⎧j) = s ⎪ π ⎪ 0≤j≤K i+j,i , ⎨ i=0 πj · e = s+K−j ⎪ ⎪ πi+j,i , K + 1 ≤ j ≤ s + K ⎩
(6)
simultaneously, the stationary probability distribution can be denoted as πk = π2s Rk−2s ,
2s < k ≤ s + K − 1
and (π0 , π1 , · · · , π2s ) is determined by the following set of equations ⎧ (π0 , π1 , · · · , π2s )B[R] = (π0 , π1 , · · · , π2s ) ⎪ ⎪ ⎪ ⎪ ⎨ 2s−1 πi e + π2s (I − RK−s )(I − R)−1 + πs+K e = 1 i=0 ⎪ ⎪ = (πs+K−1 Ps+K−1,s+K + πs+K−2 Ps+K−2,s+K ) π ⎪ ⎪ ⎩ s+K (I − Ps+K,s+K )−1 (7) The proofs of Theorem 4.2 and Theorem 4.3 are similar to literature [6,9]. Theorem 4.2 and Theorem 4.3 indicate that K affects the structure of stochastic matrix B[R], hence, also affects the value of stationary probability vector. However, it is hard to obtain analytic expression of rate matrix R from matrix equation (2), Latouche[7] developed valid recursive method to resolve R: R0 = 0, Rn+1 = −
s+2 i=2
−1 Rni P2s,2s+2−i + P2s,2s+2 P2s,2s+1 .
(8) then Rn converges to the smallest nonnegative solution R of matrix Eq.(2).
5 Analysis of performance indices (1) Under the equilibrium condition, the stationary probability distribution of the first class of customers in the system are as follows
(9)
i=0
The mean of the first class of customers in the system is s+K ¯ (1) = j · P (L(1) = j) = L s+K−j j=1 k s s+K j j πi+j,i + πi+j,i j=1
i=0
(10)
i=0
j=K+1
On the basis of the characteristic of discrete-time Erlang loss system, we know that the conditions which the first class of customers have rejected in the instant t = n− are as follows: a new customer of the first class arrived when there are s + K customers in the system, at the same time, there is no customer departure in this slot. Therefore, the loss probability of the first class of customers is (1)
PLoss =
s
µ ¯s−i ¯i2 πs+K,i 1 µ
(11)
i=0
(2) Under the equilibrium condition, the stationary probability distribution of the second class of customers in the system is as follows P (L(2) = j) =
s+K
πi,j ,
0≤j≤s
(12)
i=0
The mean number of the second class of customers in the system is
¯ (2) = L
s j=1
j · P (L(2) = j) =
s j=1
j·
s+K i=0
πi,j
(13)
MA Zhan-you, et al./Systems Engineering – Theory & Practice, 2007, 27(1): 91–98
Similarly, the loss probability of the second class of customers is (2)
PLoss = µ ¯s2 ·
s+K
πi,s
(14)
i=s
(3) The mean utilization ratio of servers is ¯ (2) , s}/s ¯ (1) + L E = min{L
(15)
6 Numerical examples Assume that there are two classes of customers, one is the first class (with delayed and loss) and the other is the second class (with loss). The arrival rate of the second classes of customers changes from 0 to 0.4 and the service rate of the second classes of customers changes from 0.01 to 0.4, the ratio of arrival rate and service rate of two classes of customers are both 1:5 (most of customers are with loss). Now, we analyze the performance indices of queue with two classes of customers when the numbers of servers are 4, 8, 16 and 32. Various performance indices can be seen in Table 1 when the loss probability of the first class of customers in the system equals 10% (the service rate of the second class of customers is 0.04). Figure 2 shows how the number of the first class customers changes with the arrival rate of the customers. It can be observed that when the arrival rate increases, i.e. the system offered load increases, the number of the first class customers increases sharply. When the arrival rate is large (the arrival is fixed), the increased number of servers will considerably decrease the number of the customers for the first class. Figure 3 shows the number of the second class of customers as a function of the service rate. From Figure 3, it can
be seen that when the service rate increases (higher the service capability), the number of the second class of customers decreases; simultaneously, when the service rate is low, the change for the number of the servers will change the number of the second class customers. Figure 4 shows the relationship between the loss probability for the first class of customers and the arrival rate. The loss probability for customers increases with the system offered load. On the other hand, the increased number of the servers will considerably decrease the loss probability of the customers. In Figure 5, we plot the loss probability for the second class customers as a function of the service rate. It can be noted that both the increased number of the servers and increased service capability can have the effect of decreasing the loss probability of the system. In Figure 6, we plot the mean number of the first class of customers as a function of the arrival rate and waiting space K (where s=8). It can be noted that when the system offered load increases, the increased number of the servers considerably decreases the number of the first class of customers. In Figure 7, we plot the mean number of the first class of customers as a function of the service rate of the second class customers and waiting space K (where s=8). It can be noted that the effect of increasing the service capability is larger for the number of the first class customers than the effect of increasing the waiting space. Figure 8 and Figure 9 show the relationship between the utility probability of servers and the arrival rate of the second class of customers, and between the utility probability of servers and the service rate respectively. It can be noted that both the increased system offered load and the increased service capability can have the effect of the utility probability of servers.
Table 1. The Loss Probability for the First Class of Customers in the System Equals 10%
Number of servers
4 8 16 32
Waiting Loss probability space of the first class K of customers 10 10 10 10
Loss probability of the second class of customers
10% 10% 10% 10%
0 0 0.1% 0.3%
Figure 2. Number of the first class of customers vs. arrival rate of the second class of customers
Ratio of arrival rate of two classes of customers are 1:5 0.02 0.08 0.16 0.32
Number of the first class of customers
Number of the first class of customers
0.1268 0.4109 1.0150 2.4011
1.8156 5.3034 12.5200 27.4940
Mean utilization ratio of servers 48.56% 71.48% 84.59% 93.41%
Figure 3. Number of the second class of customers vs. service rate of the second class of customers
MA Zhan-you, et al./Systems Engineering – Theory & Practice, 2007, 27(1): 91–98
Figure 4. Loss probability of the first class of customers vs. arrival rate of the second class of customers
Figure 8. Utility probability of servers vs. arrival rate of the second class of customers
Figure 5. Loss probability of the second class of customers vs. service rate of the second class of customers Figure 9. Utility probability of servers vs. service rate of the second class of customers
7
Figure 6. Number of the first class of customers vs. the wait space and arrival rate
Figure 7. Loss probability of the first class of customers vs. wait space and arrival rate
Conclusions
Priority queues have important application background. In this study, we analyzed a class of queues with preemptive resume priority, and obtained the stationary probability distribution of joint queue length of two classes of customers. Moreover, we obtained mean queue length of all the classes of customers, and obtained loss probability and utility ratio of the system. From these indices, we can improve the design of capacity of waiting space to satisfy requirements of different customers.
References [1] Douglas R M. Computation of steady-state probability of M/M/1 priority queues. Operation Research, 1981, 29(5): 945–948. [2] Andreas B, Manfred B. On a two-queue priority system with impatient and its application to a call center. Methodology and Computation in Applied Probability, 1999, 1: 191–210. [3] Hu G. Performance analysis of communication network on the M/M/s queue with two priorities. Jiangsu University transaction (Sciences), 2002, 23(4): 91–94. → [4] Niu Z, Kawai T, Akimaru H. Analysis of PH-MRP− + M/M2 , M1 /s(0, ∞) mixed loss and delay system with partial preemptive priority. IEICE Technical Report, 1990: 81–90. [5] Sun R, Li J. Foundation of queueing system. Beijing: Sciences Publishing house, 2002. [6] Tian N. Stochastic service system with vacation. Beijing: Beijing University Publishing house, 2001.
MA Zhan-you, et al./Systems Engineering – Theory & Practice, 2007, 27(1): 91–98
[7] Latouche G, Ramaswami V. Introduction to matrix analytic methods in stochastic modeling. Philadelphia: Society for Industrial and Applied Mathematics, 1999. [8] Neuts M. Matrix-geometric solutions in stochastic models. Baltimore: Johns Hopkins University Press, 1981.
[9] Tian N, Yue D. Quasi birth-and-death process and matrixgeometric solution. Beijing: Sciences Publishing house, 2002. [10] Hunter J. Mathematical techniques of applied probability. Academic Press, 1983, 2.