CC BY
36
Single Server Retrial Queueing System with Catastrophe
Neelam Singla 1 and Ankita Garg2* •
1 Department of Statistics, Punjabi University, Patiala-India, [email protected] 2* Department of Statistics, Punjabi University, Patiala-India, [email protected]
Abstract
The present paper analyses a retrial queueing system with Catastrophe. Primary and secondary customers follow Poisson processes. Inter arrival and service times are Exponentially distributed. Catastrophe occurs on a busy server and follows Poisson process. The server is sent for repair after its failure. The repair times are also Exponentially distributed. Steady state and time dependent solutions for number of customers in the system when the server is idle or busy are obtained. The probability of the server being under repair is obtained. Some performance measures are also evaluated. Numerical results are obtained and represented graphically.
Keywords: Queueing, Retrial, Catastr ophe, Repair .
1. Introduction
We have seen in many real life situations that sometimes a customer on arrival does not get the service instantly. So he tries for the service after some random amount of time which is popularly known as retrial. Retrial queue is a model of this kind of system if the server is not free, the customer leaves the service area and joins the virtual queue known as orbit. Thereafter it retries from the orbit after a random amount of time to get service. The queueing systems with these repeated attempts have been used in many field such as telecommunication, computer netw orks, data transmission, etc. The analysis of such systems lead to the identificatio of a new class of queueing systems known as retrial queueing systems.
For example: In call centers wher e when customers call, if they are able to reach a live negotiator immediately , they are answered else they repeat the call after a couple of minutes. The work on retrial queues in its early stages can be found in [1]. In [2] the author discussed some important single server retrial queueing models and represented analytic results. In [3] the single server retrial queue with finit number of sources is analyzed and customer 's arrival distribution, busy period and waiting time process is established. Time dependent probabilities for exact number of arrivals and departur es from the system when the server is free or busy are obtained in [4]. An explanation of the retrial queueing system is shown in the following diagram.
Secondary Arrivals
Figure 1: Basic Structure of Retrial Queueing System
Recently, a new concept of catastrophe has been introduced in queueing systems. The word catastr ophe refers to a sudden, unexpected failur e of a machine, computer netw ork, electr onic system, communication system, etc. Catastr ophes occur randomly, eradicating all customers present in the system and temporarily inactivating the service facilities. Catastrophe resets the system from current state to zero state at random time intervals. Catastrophe may come from outside the system or from another service station. Retrial queueing models with catastr ophe have applications in call centers, computer netw orks and in telecommunication systems that depend on satellites. In population dynamics, catastr ophe can be consider ed as the natural disasters such as floods storms, etc. On the other hand when we talk of catastrophe in queueing systems, it deletes all the customers present there and causes breakdo wn of the server. A basic example of retrial queueing system with catastr ophe is in call centers wher e if customers are able to reach a live negotiator immediately upon making a call, they are answered else they repeat the call after a couple of minutes. Furthermore, loss of all the customers and inactiv ation of the server will take place as a result of an incidental power failur e or a virus attack. The diagram below shows the retrial queue with catastr ophe.
Secondary Arrivals
Figure 2: Basic Structure of a Retrial Queueing System with Catastrophe
The initial work on catastrophe occurring in a simple Markovian queue could be referred from [5] and [6]. In [7] the author discussed mean queue length and the asymptotic behavior of the probability of server being free. Also the steady state probabilities are obtained. The transient solution for the system with ser ver failur e and non-zer o repair time on M/M/1 queueing system with catastr ophe is obtained by [8] .
In this research paper when a server fails, it is sent for repair immediately and after getting repaired, the server comes back to its working position and the system becomes ready to accept new customers.
The paper has been organized in the following sections.
In section 2 the complete mathematical description of the model is defined Also, the difference-differential equations are derived in this section. Steady-state solution of the model along with the expected number of customers in the system is given in section 3. In section 4 the transient state probabilities and the probability of server being under repair are obtained. In section 5 verificatio of results is given. The numerical results are obtained and represented graphically in section 6. In section 7 the busy period probabilities of system and the server are obtained numerically and presented graphically . Section 8 discusses the conclusion and finall the references are listed.
2. Model Description
In this paper , a single ser ver retrial queueing system with catastr ophe is consider ed. In this system, the customers arrive according to Poisson process. On arrival if a customer find the server busy, he joins the orbit and retries from the orbit. These retrials are consider ed to be secondar y arrivals. Catastrophe occurs on a busy server following Poisson process. It is assumed that the catastrophe occurs only when the system is non-empty and the server is busy. It has no effect on the system when the system is empty . Catastr ophe makes system empty and also causes breakdown of the server. Once the system becomes empty and the server breaks down, it
is sent for repair immediately . Further, it is assumed that during the repair time no arriv al can take place.
Assumptions : The assumptions underlying the model are listed below.
• Arrival Process: The primary customers arrive at the system according to Poisson process with mean arrival rate A.
• Retrial Process: The secondar y customers arrive at the system according to Poisson process with mean arrival rate 9.
• Service Process: The service times are Exponentially distributed with parameter p.
• Catastr ophe: Catastr ophe occur at the system accor ding to Poisson pr ocess with rate £.
• Repair: The repair time is Exponentially distributed with parameter t.
The input fl w of primary calls, intervals between repetitions, service times, catastrophes and repair times are statistically independent.
Shift operator E is used to increase the value of argument x by h so that Ef (x) = f (x + h), E2f (x) = E[Ef (x)] = E[f (x + h)] = f (x + 2h) and so on. Here h is the equal interval of spacing. Laplace transfor mation f(s) of f (t) is given by:
c TO
-st,
f (s) = e- f (t)dt, Re(s) > 0;
Jo
The Laplace inverse of
Q( p)
tmk-leakt d1-1 (Q(p)\ . m w , r ■ ,1
P(p)- E E m - /);(/ - i)ix dPM{p(P^)(p -ak)k Vp =flfc, fli =akfori =k
wher e
P(p) = (p - ai )mi (p - «2)m2.......(p - an)mn
Q(p) is a polynomial of degr ee < m1 + m2 + m3 +............ + mn - 1
If L-1 {f(s)} = F(t) and L-1 {g(s)} = G(t) then
L-1 {f (s)g(s)} = t F(u)G(t - u)du = F * G Jo
F * G is called the convolution of F and G.
2.1. Notations
Pn,0(t) = Probability that there are n customers in the system at time t and the server is free. Pn,1 (t) = Probability that there are n customers in the system at time t and the server is busy. Q(t) = Probability that the server is under repair at time t. Pn(t) = Probability that there are n customers in the system at time t.
Pn(t) = Pn,0 (t)+ Pn,1 (t) V n > 0; and Pn,1 (t) = 0 for n = 0;
Initially
P0,0 (0) = 1; Pn,0 (0) = 0, n = 0; P^ (0) = 0, Vn; Q(0) = 0;
2.2. The Difference-Differential equations governing the system are:
|Ao (t) = -(A + nd )Pn,0 (t) + pPn+1,1 (t) dtAo (t) = -APo,o (t) + TQ(t) + (t)
n > 1 (1) (2)
dtPn,i (t) = -(A + V + £)Pn,i (t) + APn-1,0 (t)+ APn-i,i (t)(l - Sn>i ) + ndPn,o(t) n > 1 (3)
j-tQ(t) = -TQ(t) + £ E Pn,i(t)
n=i
wher I
^n,i =
(4)
i, when n = i 0, other wise
3. The steady-state difference equations governing the system
Taking Pn(t) ^ Pn and ftPn(t)^ 0 as t ^ œ
(A + nd)Pn,o = V-Pn+ii
APo,o = vPi,i + t Q (A + V + £)Pn,i = APn-i,o + APn-i,i (i - Sn,i ) + ndPn,o
œ
tQ = £ E Pn,i
n=i
wher e
^n,i =
i, when n = i o, other wise
n>i ni
(5)
(6)
(7)
(8)
3.1. Steady-state solution of the problem
Using Ef(x) = f (x + 1), equations (5) and (7) are represented as
[(A +(n + i)d)E]Pnfi - VE2Pn,i = o [A + ((n + i)d)E]Pnfi + [A - (A + V + £)E]Pn,i = o
In order to fin the solution of the above system of equations, we need
n>2 n2
(9) (io)
E[v(n + i)dE2 - (A(A + (n + i)d + £) + d((n + i)(v + £)))E + A(A + (n + i)d)] = o n > 2
(ii)
The values of Pn 0 and Pn1 are given by
2
n
izi
Pn,o = E aiz
i=o
2
Pn,i = E b
n
izi
n2
n2
i=o
œ
wher e z0,z1, z2 are the roots of (11) with z0 = 0 and ai,bi, i = 0,1,2 are the arbitrar y constants to be evaluated. Other two roots of the equation (11) are given by
Z1, Z2
2p(9+n9)
A(A + (n + 1)9 + £) + 9((n + 1)(p + £)) ±
292 (p£ - pA + £A)
+ 2n9 (p2 + £ + A2 + 2p£ - 2pA + 2£A) + 2A2 (29£ - p9) + 2nA2 (29£ + 9A - 9p)
+ 29p£A(n + 1) + A2 (£ + 9 + A2) + 2A3 (9 + £) + 2£2 9A(n + 1) + 92 (p2 + £2)
1/2
(12)
Clearly the root z1 is always greater than 1 and the root z2 is always less than 1. For the convergence of a solution, a root greater than or equal to 1 must be rejected. So when Zi > 1, ai and bi are taken as equal to zero. Since here z1 is greater than 1, so we take a1 = b1= 0
As Z0 = 0 and a1 = b1= 0, therefore the values of Pn0 and Pn1 are given by
Pn,0 = a2n > 2 Pn,1 = b2zn n > 2
From equations (5) and (7) for n=1 the probabilities P10 and P11 are given by
P1,0 = TT^T P2,1
p
P1,1
(A + 9) A
(A + 9)
P0,0 +
b2z22
p9
-Jb2 z\
A + p +10,0 ' (A + p + £)(A + 9)
By substituting the above values in equation (7) for n = 2 and for n = 3, we have pA , 9 V 2 . A2 29
b2z22 =
(A + 9)(A + p + £ )\ A + p + £
b2 z2
A
+ ■
A + p + £ 39
1 +
29
a2 z2 +
pA2
A + p + ' (A + 9)(A + p + £)2
3
1 +
9
A + p + £
b2z22
-a2z2 +
A
A + p + £ 2 \A + p + £ On solving equations (17) and (18) we get
P
0,0
a2
A
A
B-A
2AC - B (xrpT£ (1 + 2C) + 3z2c)
P
1
2C
A
A
A
0,0
+
A
b2 = -<
Bz2l2AC - B(A+p^ (1 + 2C)+ 3z2c)
and the value of P0,0 can be found by using the relation
TO
P0,0 + E (Pn,0 + Pn,1 ) + Q = 1
n=1
A + p + £
P
0,0
(13)
(14)
(15)
(16)
9 2 A2 29 2
1 + m—m 2 z2 + Tin—P00 + m—m a2 z2 (17)
v A + p + £J 2 (A + p + £ )2 A + p + £ 2
(18)
(19)
(20)
1
2
z
2
B
2
After simplificatio
P
0,0
1 + ( 1 + T
A
+
A + V + £J \A + 0
V
+
1 + T
T
V0
2C
A
A+V+t
B
A
A+V+t
A
1
+ ^
2ABC - B2 (1 + 2C) + 3z2c) B Va + v +t
(A + 0){A + V + t)
2
+
+
1 - z2
1 + T
1 - z2
A
A+V+t
A
A+V+t
BA
{2AC - b(■A+A+t (1 + 2C)+ 3z2C))
2C
a
A+V+t
A
A+V+t
A
1
+ TT
2 ABC - B2( (1 + 2C)+ 3z2 c) B Va + V + t
A
2
-1
(21)
wher e
A = z2 -B=1-C
VA2
(A + 0)(A + V + t )2 VA
1+ 1+
0
A+V+t 0
(A + 0)(A + V +1 )\ A + V + t 0
A+V+t
Hence by using the values of a2, b2 and P0 0, the probabilities Pn 0, Pn1 and Q are completely known for various values of n.
3.2. Expected number of customers in the system
Expected number of customers in the system is given by
Ls = Ls,0 + Ls,1
wher e
Ls,0 denotes the expected number of customers in the system when the server is free. Ther efor e, by definitio of expectation
TO
Ls,0 = E npn,0
n = 1
TO
= P1,0 + E nPn,0
n=2
P1,0 + a2 E
nz
n=2
= P1,0 + a2 z2
L(1 - z2 )
1
(22)
Similarly
Ls,1 denotes the expected number of customers in the system when the server is busy.
Ls,1 = P1,1 + b2 z2
L(1 - z2)2
1
(23)
2
A
2
1
2
B
1
1
By using (24) and (25)
Ls = Pi,0 + Pi,i + (a2 + b2 = Pj + (a2 + b2 )
(1 - z2 )2
- 1 z2
(1 - z2 )2
- 1 z2
(24)
4. Laplace transform of difference-differential equations
Using the Laplace transfor mation f(s) of f (t) given by
f(s)
e f (t)dt, Re(s) > 0;
in the equations (1)-(4) along with the initial conditions, we have
wher e
n1
(s + A + n6)Pn,0 (s) = pPn+1,1 (s) (s + A)Po,o(s) - 1 = tQ(s) + (s) (s + A + v +£)Pn,1 (s) = APn-1,0 (s) + APn-1,1 (s)(1 - Sn>1 ) + n9Pnfi (s) n > 1
L n-œ
(s + T)Q(s) = £ E PnA (s)
n=1
$n,1 =
1, when n = 1 0, other wise
(25)
(26)
(27)
(28)
4.1. Transient solution of the Problem
Solving equations (25)-(28) recursiv ely, we have
Po,o (s)
Pn,0 (s)
1 Q (s) + 7-^T P1,1 (s)
(s + A) (s + A)
(s + A)
s + A + n0
A
n+1
y ,
Lkti V s + A + v + £
n-k+1
wher i
1
Vk (s)Pk,o (s) +
if k = 1
A
s + A + v + £ n1
Pu (s)
(29)
(30)
Vk(s)=\ 1 + Ï+A+ÏH if k = 2 ton
.ÏTÂ+ÏH if k = n + 1
Pu (s)
A
s + A + V + £
1 Q (s) + T^ Pu (s)
Pn,1 (s) = y k=1
wher e
A
s + A + v + £
(s + A) (s + A)
n-k
Vk (s)Pk,0(s)
(s + A)
+
0
Vk(s)
+
A
s + A + v + £
s + A + v + £
n-1
Pu (s)
n2
P1,0 (s) (31)
(32)
1 +
k0
if k = 1
if k = 2 to n 1
s+A+v+£
k0 if k = n
s+A+v+£
1
œ
n
n
1
X TO
Q{s) = (sh P« (»)
Taking the Inverse Laplace transfor m of equations of (29)-(33), we have
Po,o (t) = e—M + te—xt * Q(t) + pe—xt * Pu (t)
(33)
(34)
Pn,0 (t) = pXne—(x+ne)t
1 n—1 lT 1
J__e — (p+Z)t Y i__1
L(p + Z )n
[(p + Z )n—k+ 1
1 n—k
1 — e—(p+z)t y
t=0 t! (p + z)n—t\ 1
k=2
t! (p + Z)
t=0
n—k+1 tT
n—k—r+1
__e—(p+z)t y
[(p + Z)n—k+2 y t! (p + Z)n—k—T+2\
- 1 e—(p+Z)t -
.(p + Z) — (p + Z). * P1,1 (t)
* P1,o (t) + e—(x+ne)t y pAn—k+1
* Pk,0 (t) + e—(x+ne)t (x+ne
E(pke)x
k=2
n—k+1
* Pn+1,0 (t) + pXne—(x+ne)t
* Pk,o(t) + e—(x+ne)t(n + 1)pe
1 n — 1 .t 1
i__e—(p+Z£__1
L(p + Z )n
e0 t! (p + z)n—t\
n1
(35)
in—2
n1
Pn,1 (t) = xn—1 e—(x+P+Z)t _ * P1,0 (t) + e—(x+F+Z)f y An—k(
t
n—k—1
(n — 2)! 1,0 ^ 1 k=2 (n — k — 1)!
n —1 tn—k tn—2
y (ke)xn—k(n—y^ * Pk,0 (t) + nee—(A+p+Z)f * Pn,0 (t) + xn—1 e—(A+p+Z)f * P1,1 (t)
* Pk,0 (t) + e—(A+P+Z)t
n2
P1,1 (t) = Xe—(x+P+Z)t * Po,o (t) + ee—(x+P+Z)t * P1,0 (t)
TO
Q(t) = Ze—Tt E Pn,1 (t)
n=1
(36)
(37)
(38)
5. Verification of Results
Summing equations (29)-(33) over n we get,
TO 1
y [4,0 (s) + P„,i (s)] + Q(s) = -
n=0 s
and hence
y [Pn,0 (t) + Pn,1 (t)] + Q(t) = 1
n=0
which is a verificatio for our results.
6. N umerical Solution and Graphical Representation
The Numerical results are generated using MATLAB programming for the case p = (x) = 0.8
y = (p) = 0.9, t = (p) = 0.6, Z = (p) = 0.4. In following tables, we list some significan probabilities at various time instants.
Table 1: At time t = 1
t P0,0 P1,0 P2,0 P3,0 P4,0 P5,0 P1,1 P2,1 P3,1
1 0.5904 0.0216 0.0033 0.0004 0.0001 0 0.226 0.0701 0.0161
TO
P4,1 P5,1 Sum
0.0029 0.0005 0.9314
Table 2: At time t = 5
t P0,0 Pl,0 P2,0 P3,0 P4,0 P5,0 P1,1 P3,1 P4,1
5 0.3644 0.0556 0.0187 0.0069 0.0032 0 0.1572 0.0489 0.0247
P5,1 Q(t) Sum
0.0145 0.2126 0.9067
Table 3: At time t = 15
t P0,0 P1,0 P3,0 P4,0 P5,0 P1,1 P2,1 P3,1 P4,1
15 0.3571 0.0528 0.0078 0.0041 0 0.1514 0.0897 0.0491 0.0274
P5,1 Q(t) Sum
0.018 0.2237 0.9811
Table 4: At time t = 25
t P0,0 P1,0 P2,0 P3,0 P4,0 P5,0 P1,1 P2,1 P3,1
25 0.3571 0.0528 0.0189 0.0078 0.0041 0 0.1514 0.0897 0.0491
P4,1 Q(t) Sum
0.0274 0.2237 0.982
Table 5: At time t = 40
t P0,0 P1,0 P2,0 P3,0 P5,0 P1,1 P2,1 P3,1 P4,1
40 0.3571 0.0528 0.0189 0.0078 0 0.1514 0.0897 0.0491 0.0274
P5,1 Q(t) Sum
0.018 0.2237 0.9959
In the following figu es, the probabilities are graphed against time. Figure 3 to figu e 5 are plotted for the case p = 0.8, q = 0.9, t = 0.6, £ = 0.4.
Figure 3: Probabilities P0,0 and Q(t) against average service times t
In figu e 3, the probabilities P00 and Q(t) (probability of the server being under repair) are plotted against time t. From the graph, we observe that the probability P00 decreases rapidly from its initial value 1 at t = 0 and thereafter becomes steady. On the other hand, the probability Q(t) increases in the beginning and then becomes stable.
Figure 4: Probabilities P10, P2,0 and P3,0 against average service times t
Figure 5: Probabilities P1,1, P2,i, P3,1 and P4,1 against average service times t
The probabilities P10, P2,0, P3,0 are plotted against time t in figu e 4 and the probabilities P1,1, P2,1, P3,1, P4,1 are plotted against time t in figu e 5. The graphs clearly indicate that all the
probabilities increase rapidly in the beginning, gradually decline to a certain extent and finall become stable. Also, it can be seen from figu e 4 that the probability of more customers in the system achieves a lower highest value. In addition, figu e 5 shows that the probability of server being busy attains a lower steady value when there are more customers in the system.
Figure 6: Effect of change in on probability Q(t)
In figu e 6, we study the effect of change in £ (catastrophe rate per unit service time) on the probability Q(t)(probability of server being under repair). From the graph it can be seen that whenever the catastrophe rate per unit service time increases, the probability Q(t) also increases which is as desir ed.
t(aver'age service tuiins)—>
Figure 7: Effect of change in t on probability Q(t)
In figu e 7, the effect of change in t (repair rate per unit service time) on the probability Q(t) is studied. From the graph it is clearly visible that whene ver the repair rate per unit service time increases, the probability Q(t) decreases.
7. Busy Period Probabilities
This section discusses the busy period probabilities of server and system. In terms of probability , busy servers are deter mined as follows:
P(Serveris busy) = y Pn,1 (t) (39)
n>1
And busy systems are determined as follows:
P(System is busy )= y (Pnfi(t) + Pn,i (t)) + Q(t) (40)
n>0
7.1. Numerical and Graphical Representation of Busy Period Probabilities
The numerical results are obtained using MATLAB programming and following [9]. The Probabilities of system busy and server busy are obtained for different values of p keeping the other parameters constant and are presented in the table below:
Table 6: Probabilities of system busy and server busy for different values of p
t Probability(System Busy) Probability(S erver Busy)
p = 0.4 p = 0.6 p = 0.8 p = 0.4 p = 0.6 p = 0.8
0 0 0 0 0 0 0
1 0.1916 0.2711 0.341 0.1833 0.2548 0.3156
2 0.2361 0.3283 0.4042 0.212 0.2853 0.3434
3 0.2514 0.3463 0.4203 0.2169 0.2884 0.3432
4 0.2577 0.3525 0.4235 0.2178 0.288 0.3407
5 0.2604 0.3544 0.4229 0.2179 0.2872 0.3387
6 0.2615 0.3548 0.4217 0.2177 0.2865 0.3373
7 0.262 0.3547 0.4207 0.2176 0.2861 0.3366
8 0.2621 0.3546 0.4201 0.2175 0.2859 0.3361
9 0.2622 0.3545 0.4197 0.2175 0.2857 0.3359
10 0.2622 0.3544 0.4194 0.2174 0.2857 0.3357
A graph depicting the probabilities of server and system busy is also included.
Figure 8: Probabilities of system busy and server busy against average service times t
In figu e 8, the probabilities of system busy and server busy are plotted against time t for the case p = 0.8, n = 0.9, t = 0.6, £ = 0.4. It is clearly visible from the graph that probability of system busy is higher than probability of server busy. Both probabilities increase rapidly in the beginning and then become stable with time.
Figure 9: Probability of system busy for different values of g
Figure 10: Probability of server busy for different values of g
In figu es 9 and 10, the probability of system busy and the probability of server busy are plotted respectively against time t for different values of g keeping other parameters constant. As we know when catastrophe occurs, the server breaks down and system becomes empty. So the probability that both system and server remain busy attains lower value for greater values of catastr ophe rate per unit ser vice time.
8. Conclusion
We have modeled a single server retrial queueing system with catastr ophe to quantify various perfor mance measur es and understand characteristics of related systems. The catastr ophe has significan impact on businesses, computer netw orks, etc. It is very important to manage the risk of catastr ophe for the smooth functioning of the system. In this paper, the steady-state and transient state probabilities for the number of customers in the system when the server is busy or idle are obtained by solving difference-differential equations. In addition, the probability that server is under repair is also obtained. Some performance measur es are given. Numerical solutions and busy period probabilities are obtained by using MATLAB programming and presented graphically . This model is applicable in call centers, computer networks, etc.
R eferences
[1] JW Cohen. Basic problems of telephone traffic the ory and the influenc of repeated calls. Pillips Telecomm. Rev., 18(2), 1900.
[2] Tao Yang and James G. C. Templeton. A survey on retrial queues. Queueing systems, 2(3): 201-233, 1987.
[3] Gennadi I Falin and Jesus R Artalejo. A finit source retrial queue. European Journal of Operational Research, 108(2):409-424, 1998.
[4] Neelam Singla and Sonia Kalra. Performance analysis of a two-state queueing model with retrials. Journal ofRajasthan Academy of Physical Sciences, 17:81-100, 2018.
[5] Xiuli Chao. A queueing network model with catastrophes and product form solution. Operations Research Letters, 18(2):75-79, 1995.
[6] Antonio Di Crescenzo, Virginia Giorno, Amelia G Nobile, and Luigi M Ricciar di. On the m/m/1 queue with catastrophes and its continuous approximation. Queueing Systems, 43(4): 329-347, 2003.
[7] B Krishna Kumar and D Arivudainambi. Transient solution of an m/m/1 queue with catastr ophes. Computers & Mathematics with applications, 40(10-11):1233-1240, 2000.
[8] B Krishna Kumar , Achyutha Krishnamoorthy , S Pavai Madhesw ari, and S Sadiq Basha. Transient analysis of a single ser ver queue with catastr ophes, failur es and repairs. Queueing systems, 56(3):133-141, 2007.
[9] Brian D Bunday. Basic queueing theory. Arnold, 1986.
