An enhanced algorithm to solve multiserver retrial queueing systems with impatient customers
Full text
An Enhanced Algorithm to Solve Multiserver Retrial Queueing Systems with Impatient Customers Tien Van Do a,∗Nam H. Do bJie Zhang c aMTA-BME Information Systems Research Group, Department of Networked Systems and Services, Budapest University of Technology and Economics, H-1117, Magyar tud´osok k¨or´utja 2., Budapest, Hungary. bInter-University Centre for Telecommunications and Informatics, Budapest University of Technology and Economics, 4028 Debrecen, Kassai ´ut 26., Hungary cCommunications Group, the Department of Electronic and Electrical Engineering, University of Sheffield, Mappin Street, Sheffield, S1 3JD UK Abstract The homogenization of the state space for solving retrial queues refers to an approach where the performance of the M/M/c retrial queue with impatient customers and cservers is approximated with a retrial queue with a maximum retrial rate restricted beyond a given number of users in the orbit. As a consequence, the stationary distribution can be obtained by the matrix-geometric method, which requires the computation of the rate matrix. In this paper, we revisit an approach based on the homogenization of the state space. We provide the exact expression for the conditional mean number of customers based on the computation of the rate matrix Rwith the time complexity of O(c). We develop simplified equations for the memory-efficient implementation of the computation of the performance measures. We construct an efficient algorithm for the stationary distribution with the determination of a threshold that allows the computation of performance measures with a specific accuracy. Keywords: retrial queues, matrix-geometric method, spectral expansion, efficient algorithm T. V. Do, N. H. Do, J. Zhang. An Enhanced Algorithm to Solve Multiserver Retrial Queueing Systems with Impatient Customers. Computers & Industrial Engineering, DOI:10.1016/j.cie.2013.04.008, 2013 ∗Corresponding author Email addresses: [email protected] (Tien Van Do), [email protected] (Nam H. Accepted to Computers & Industrial Engineering
1 Introduction Retrial queues have been used to take into account a phenomenon in modern information and telecommunication systems that blocked customers may rerequest for service after a certain timeout [1–10]. In retrial queues a client who does not receive the allocation of a server joins the orbit and later initiates a request for service. The M/M/c retrial queue has been analyzed by many researchers because the stationary distribution when the number of servers is larger than two can be only obtained using approximate techniques [1,6– 8,11,12]. Falin [13] presented necessary and sufficient conditions for ergodicity of the retrial queues M/M/c. A well-known approximation is based on the truncation of the state space at a sufficiently large level related to the number of customers in the orbit [13]. Another approximation based on the homogenization of the model was pioneered by Neuts and Rao [14], where the M/M/c retrial queue is approximated by the multiserver retrial queue with the total retrial rate that does not depend on the number of clients in the orbit as long as the orbit contains the number of clients greater than the specified value N. Note that the discussion for the choice of Nis presented in the recent book by Artalejo and G´omez-Corral on retrial queues [8]. With this assumption, the stationary probabilities of the M/M/c retrial queue can be estimated by any algorithm [15–19] based on the matrix-geometric method (MGM). Recently, Domenech-Benlloch et al. [20] considered a multiserver retrial queue with the impatient phenomenon of customers waiting in the orbit. They proposed two different generalized truncated methods (called HM1 and HM2) based on the homogenization of the state space beyond a given number of users in the retrial orbit. The steady-state probabilities of the multiserver retrial queue with impatient customers are approximated with a modified retrial queue where the retrial rate beyond a certain level only depends on the conditional mean value of the number of customers in the orbit. Domenech-Benlloch et al. [20] also compared their methods with other well-known algorithms that belong to different categories [11] (approximations, finite truncated methods, generalized truncated methods). The authors [20] showed that the proposed HM2 method outperforms previous approaches from the aspect of accuracy at the price of increasing computation cost. Based on the HM2 algorithm of Domenech-Benlloch et al. [20], our contributions allow an efficient computation for the stationary distribution and the performance measures. First, we revisit an approach based on the homogenization of the state space and provide an efficient method with the time complexity of only O(c) to compute the rate matrix R. The method is based Do). 2
on a property that the characteristic matrix polynomial has only a single nonzero eigenvalue and this single non-zero eigenvalue can be computed using the bisection method. Second, we derive an exact expression for the conditional mean number of customers. Third, we develop simplified equations that allow the memory-efficient implementation of the computation of the performance measures. Fourth, we construct an efficient computation for the stationary distribution with the determination of a threshold, which guarantees a specific accuracy for the computation of performance measures. The rest of this paper is organized as follows. In Section 2 we summarize the considered queueing model with impatient customers. In Section 3 we present our new results that serve as the foundations of the computation. In Section 4 we provide some numerical results to illustrate the efficiency of our algorithm. Finally, Section 5 concludes our paper. 2 A Retrial Queueing Model with Impatient Customers We consider a retrial queueing model with chomogenous servers and impatient customers. Inter-arrival times of customers are exponentially distributed with parameter λ. Holding times are exponentially distributed with parameter µ. Random variable ג(t) represents the number of occupied servers at time t, hence 0 ≤ג(t)≤cholds. A client joins the orbit in order to wait and retry upon when ג(t) = c. Let k(t) be the number of clients in the orbit waiting for retrial at time t. Each customer retries with rate µr. Hence, the total effective retrial rate, when k(t) = j, is jµr. A retrying customer either leaves the queue with probability Pim if all servers are busy upon the retrial or rejoins the orbit with probability 1 −Pim. Note that a time between subsequent retrials of a specific user follows the exponential distribution with parameter µr. This system can be represented by two-dimensional continuous-time Markov chain (CTMC) Y={ג(t),k(t)}with state space {0,1,...,c} × {0,1,...}. Let the steady-state probabilities of CTMC Ybe denoted by πi,j = lim t→∞ Pr(ג(t) = i, k(t) = j). Define the row vector vj= [π0,j ,...,πc,j]. 2.1 Notations CTMC Yis driven by the following transitions. (a) Aj(i, k) denotes the transition rate from state (i, j) to state (k, j) (0 ≤ i, k ≤c;j= 0,1,...), which is caused by either the arrival of a customer (when i < c) or the leaving of a client after the expiry of a holding time. 3
Matrix Ajis of size (c+ 1) ×(c+ 1) with elements Aj(i, k). Since Aj is j-independent, it can be written as Aj=A. The nonzero elements of Ajare Aj(i, i −1) = iµ for i= 1,...,c+ 1, and Aj(i, i + 1) = λfor i= 0,...,c. Because Ajis j-independent, it can be written as Aj=A= 0λ0... 0 0 0 µ0λ . . . 0 0 0 . . .. . .. . .. . .. . .. . .. . . 0 0 ... (c−1)µ0λ 0 0 ... 0cµ 0 ,∀j≥0. (b) Bj(i, k) represents the one-step upward transition rate from state (i, j) to state (k, j + 1) (0 ≤i, k ≤c;j= 0,1,...), which is caused by the arrival of a request when all servers are busy (i.e., when i=c), thus increasing k(t) by 1. Matrix Bj(B, since it is j-independent) is of size (c+ 1) ×(c+ 1) with elements Bj(i, k). The only nonzero element of Bj is Bj(c, c) = λ. Thus, we get Bj=B= 000... 000 000... 000 . . .. . .. . .. . .. . .. . .. . . 0 0 ... 000 0 0 ... 0 0 λ ,∀j≥0. (c) Cj(i, k) is the transition rate from state (i, j) to state (k, j −1) (0 ≤ i, k ≤c;j= 1,2,...), which is due to the successful retrial of a request from the orbit. Matrix Cjis of size (c+ 1) ×(c+ 1) with its elements Cj(i, k). The nonzero elements of Cj(j≥1) are Cj(i, i + 1) = jµrfor i= 0,...,c and Cj(c, c) = jµrPim. Matrix Cj(∀j≥1) with elements Cj(i, k) is written as Cj= 0jµr0... 000 0 0 jµr... 000 . . .. . .. . .. . .. . .. . .. . . 0 0 ... 0 0 jµr 0 0 ... 0 0 jµrPim ,∀j≥1. 4
Note that C0= 0 by definition. Let DA,DCand DCj, j ≥1 denote diagonal matrices with the diagonal elements DA(i, i) = Pc k=0 A(i, k), DC(i, i) = Pc k=0 C(i, k) and DCj(i, i) = Pc k=0 Cj(i, k) for i= 0,...,c. The balance equations, which equate the probability fluxes from and to the states of CTMC Y, and the normalization equation pertaining to CTMC Ycan be written as follows (see [3,8]): v0Q(0) 1+v1Q(1) 2=0,(1) vj−1Q(j−1) 0+vjQ(j) 1+vj+1Q(j+1) 2=0(j≥1),(2) ∞ X j=0 vjeT= 1.0 (normalization), where Q(j) 0=B, j ≥0; Q(j) 1=A−DA−B−DCj, j ≥0; Q(j) 2=Cj, j ≥1 and eis the row vector of size c+ 1 with each element equal to unity. Using the similar argument as in [3–5,8], the infinitesimal generator matrix [15,16,21] of Y, that satisfies [v0,v1,...]QY=0, can be constructed from equations (1) and (2) as follows: QY= Q(0) 1Q(0) 00 . . . . . . . . . . . . . . . . . . Q(1) 2Q(1) 1Q(1) 00 . . . . . . . . . . . . . . . 0Q(2) 2Q(2) 1Q(2) 00 . . . . . . . . . . . . 0 0 Q(3) 2Q(3) 1Q(3) 00 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Q(j) 2Q(j) 1Q(j) 0. . . . . . . . . . . . . . . . . . Q(j+1) 2Q(j+1) 1Q(j+1) 0. . . . . . . . . . . . . . . . . . . . . Q(j+2) 2Q(j+2) 1Q(j+2) 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . It is clear that QYis a block tridiagonal matrix with •QY(j, j + 1) = Q(j) 0, j ≥0,in the upper diagonal, •QY(j, j) = Q(j) 1, j ≥0 in the main diagonal, •QY(j, j −1) = Q(j) 2, j ≥1 in the lower diagonal. 2.2 An Approximation Domenech-Benlloch et al. [20] suggested that the M/M/c retrial queue with impatient customers can be approximated by the solution of the modified multiserver retrial queue with the retrial rate µr(j) = jµrif j < N M(N)µrif j≥N , 5
where M(N) = E[J|J≥N] is the conditional mean number of customers. As a consequence, the modified multiserver retrial queue is described by a CTMC Z={גZ(t),kZ(t)}with state space {0,1,...,c} × {0,1,...}, where גZ(t) represents the number of occupied servers at time tand kZ(t) is the number of clients in the orbit waiting for retrial at time t. The steady-state probabilities of CTMC Zare denoted by e πi,j = lim t→∞ Pr(גZ(t) = i, kZ(t) = j), j≥0,0≤i≤c, and the row vectors e vj= [e π0,j,...,e πc,j], j≥0. We define the transition rate matrices associated with CTMC Zas e Aj,e A,e Bj, e B,e Cjand e Cfor j≥0. Note that we have e Aj=e A=Aand e Bj=e B=Bfor j≥0. Furthermore, e Cj=Cjfor 0 ≤j < N and e Cj=e C= 0M(N)µr0. . . 0 0 0 0 0 M(N)µr. . . 0 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 . . . 0 0 M(N)µr 0 0 . . . 0 0 M(N)µrPim , ∀j≥N. For j≥N, the balance equation of CTMC Zcan be rewritten as e vj−1e Q0+e vje Q1+e vj+1 e Q2=0(j≥N),(3) where e Q0=e B, e Q1=e A−De A−e B−De C,e Q2=e C. The coefficient matrices in the difference equations (3) are j-independent. This leads to the following solution based on the MGM (see [16]) e vj=e vN−1Rj−N+1 (j≥N−1),(4) where Ris the unique minimal nonnegative solution of the quadratic matrix equation e Q0+Re Q1+R2e Q2= 0 (see [15,16]). After the computation of R, the rate matrix, the steady-state probabilities for states 0 ≤j≤N−1 can be determined by solving the balance equations pertaining to the levels 0≤j < N and the normalization equation. 6
Algorithm 1 The HM2 algorithm M0(N) = N k= 0 repeat k=k+ 1 Compute Rmatrix based on the logarithmic reduction algorithm [16] Compute Mk(N) using equation (5) until |Mk(N)−Mk−1(N)|/Mk−1(N)< ǫM Solve for vjfor j= 0,...,N Because Rand e vNdepend on M(N), we can get the fixed-point iteration M(N) = ∞ X j=N je vje ∞ X j=Ne vje =e vN[R(I−R)−1+NI](I−R)−1e e vN(I−R)−1e,(5) where Iis the identity matrix of size (c+1)×(c+1). Hence, Domenech-Benlloch et al. proposed Algorithm 1 (called HM2) in [20]. 3 An Enhanced Algorithm The stationary distribution of CTMC Yis approximated by the steady-state probabilities of CTMC Z. Therefore, we need to compute the following quantities associated with CTMC Z: •the rate matrix R, •the conditional mean number of customers M(N), •the steady-state probabilities for states 0 ≤j≤N−1, •the estimation of N. Note that Rcan be computed by the original algorithm of the MGM [15] and further improved algorithms of MGM [16,18,19]. However, the time complexity of these algorithms is O(c3). In what follows, we provide a method to compute the rate matrix (Theorem 1) in Section 3.1. We derive the exact and simplified formula for the computation of the conditional mean number of customers in the orbit (Corollary 1). As a consequence, we can compute the rate matrix Rand the conditional mean number M(N) of customers in a very efficient way. We provide a method to 7
determine the steady-state probabilities for states 0 ≤j≤N−1 in Section 3.2. We provide the new formulae of performance measures and the relation between performance measures in Section 3.3. Next, we present our new result and our algorithm for the computation of the initial value of Nin Section 3.4. 3.1 The computation of matrix Rand M(N) In Theorem 1 we prove that the characteristic matrix polynomial has only a single non-zero eigenvalue, and it can be computed using the bisection method. As a consequence, the rate matrix Rhas a special form and a method can be constructed to compute the rate matrix Rwith the computational complexity of O(c). The property that the characteristic matrix polynomial has only a single non-zero eigenvalue allows the derivation of an exact equation for the conditional mean number of customers M(N). Theorem 1 The rate matrix Rhas all rows of elements equal to zero except the last row r= [r0, r1,...,rc], where rc=xcis the single eigenvalue of characteristic matrix polynomial Q(x, M(N)) = e Q0+e Q1x+e Q2x2in the interval (0,1) (the corresponding left-eigenvector is ψc= [ψc,0, ψc,2,...,ψc,c] with ψc,c = 1) and ri=xcψc,i for 0≤i < c. The computational complexity for rcand ψcis O(c). Proof. The steady-state probabilities of the CTMC Zare expressed as e vj= c X k=0 bkψkxj−N+1 k(j≥N−1),(6) where bkare suitable coefficients to be determined using the balance equations pertaining to rows 0 to N−1 and the normalization equation, (xk,ψk), k= 0,...,c are the left eigenvalue-eigenvector pairs of Q(x, M(N)) = e Q0+e Q1x+e Q2x2inside the unit circle. They satisfy, ψkQ(xk, M(N)) = 0;det[Q(xk, M(N))] = 0, k = 0,...,c. Since the (c+ 1) ×(c+ 1) tri-diagonal matrix Q(x, M(N)) can be expressed Q(x, M(N)) = q1,1(x)q1,2(x) 0 . . . 0 0 q2,1(x)q2,2(x)q2,3(x). . . 0 0 0q3,2(x)q3,3(x). . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 . . . qc,c−1(x)qc,c(x)qc,c+1(x) 0 0 . . . 0qc+1,c(x)qc+1,c+1(x) , 8
where q1,1(x) = −(λ+M(N)µr)x, qi,i(x) = −(λ+M(N)µr+ (i−1)µ)x (i= 2,...,c), qc+1,c+1(x) = λ−(λ+cµ +M(N)µrPim)x +M(N)µrPimx2, qi,i+1(x) = λx +M(N)µrx2(i= 1,...,c), qi+1,i(x) = µix (i= 1,...,c). It is easy to verify that Q(x, M(N)) has czero-eigenvalues. Let the nulleigenvalues be x0, . . . , xc−1with corresponding independent left-eigenvectors ψ0= [1,0,...,0], ψ2= [0,1,0,...,0],. . . ,ψc−1= [0,0,...,1,0], respectively. As a consequence, Q(x, M(N)) should have a single non-zero eigenvalue xc strictly inside the unit disk to ensure that the stationary distribution of CTMC e Yexists. Let L(x, M(N)) and U(x, M(N)) denote the component matrices in the LU decomposition of Q(x, M(N)) = L(x, M(N))U(x, M(N)) for any specific value x. Due to the tri-diagonal structure, the component matrices of the LU decomposition of Q(x, M(N)) can be written as follows L(x, M(N)) = l1(x, M(N)) 0 0 ... 0 0 µx l2(x, M(N)) 0 ... 0 0 . . . . . . . . . . . . . . . . . . 0 0 . . . µ(c−1)x lc(x, M(N)) 0 0 0 . . . 0µcx lc+1(x, M(N)) , U(x, M(N)) = 1u1(x, M(N)) . . . 0 0 0 0 0 1 u2(x, M(N)) . . . 0 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 . . . 0 1 uc(x, M(N)) 0 0 . . . 0 0 1 . By equating the corresponding elements of Q(x, M(N)) and L(x, M(N)) · U(x, M(N)), and using some algebraic simplifications, we get 9
Proposition 1 The nonservice probability is expressed as Pns =λ−1Pimµra2=λ−1Pimµr N−2 X m=1 me πc,m +(M(N)−1)e πc,N−1 1−rc!.(26) Proof. Substituting (16) to the definition of a, we obtain a=e vNR(I−R)−1+NI(I−R)−1 =e vNR(I+R 1−rc ) + NI(I−R)−1 =e vNR 1−rc +NI(I−R)−1(using (13)) =e πc,N−1rR 1−rc +NI(I−R)−1(using (14)) =e πc,N−1rcr 1−rc +Nr(I−R)−1(using (13)) =e πc,N−1M(N)r(I−R)−1(using (12)) =M(N)e vN(I−R)−1(using (13)) (27) =e πc,N−1M(N)rI+R 1−rc(using (16)) =e πc,N−1M(N)r+rcr 1−rc(using (14)) =e πc,N−1M(N) 1−rc r.(28) Substituting (28) into (25), we obtain 16
a2= N−1 X m=0 me vmz+a z = N−1 X m=0 me vmz+e πc,N−1M(N) 1−rc r z = N−1 X m=1 me πc,m +e πc,N−1M(N)rc 1−rc = N−2 X m=1 me πc,m +e πc,N−1"N−1 + rcM(N) 1−rc# = N−2 X m=1 me πc,m +(M(N)−1)e πc,N−1 1−rc (using (12)). (29) Equation (29) yields (26).✷ Proposition 2 The mean number of users in the retrial orbit is Nret =a1+a2.(30) Proof. From the definition of a1and a2, we get a1+a2= N−1 X m=0 me vm(o+z) + M(N)e vN(I−R)−1(o+z) = N−1 X m=0 me vme+M(N)e vN(I−R)−1e. Utilizing (27) and the definition of Nret, we obtain equation (30). ✷. Proposition 3 We can obtain the blocking probability Pbas follows: Pb= N−2 X m=0 e πc,m +e πc,N−1 1−rc .(31) Proof. Utilizing (16), we obtain 17
Pb= N−1 X m=0 e vmz+e vN(I−R)−1z = N−1 X m=0 e vmz+e vN−1R(I+R 1−rc )z = N−1 X m=0 e vmz+e vN−1 R 1−rc z = N−1 X m=0 e vmz+e πc,N−1 r 1−rc z = N−2 X m=0 e vmz+e vN−1z+e πc,N−1 r 1−rc z = N−2 X m=0 e πc,m +e πc,N−1+e πc,N−1 rc 1−rc = N−2 X m=0 e πc,m +e πc,N−1 1−rc .✷. Proposition 4 The following relation exists between the performance measures Nret =λ µr Pns(1 −Pim) Pim +Pb!.(32) Proof. From Pds =λ−1µra1,Pns =λ−1Pimµra2,Pb=Pds +Pns and Nret = a1+a2, we get Nret =λ µr (Pns Pim +Pds) =λ µr (Pns Pim +Pb−Pns) =λ µr Pns(1 −Pim) Pim +Pb!.✷ 3.4 An Estimation of N Domenech-Benlloch et al. [20] presented some numerical results concerning choosing the appropriate value of threshold Nto achieve the required accuracy of the approximation of performance measures. However, the authors [20] did not present a systematic way to find the appropriate value of threshold N. In this section, we will show an efficient method to estimate threshold N. 18
0.0 0.2 0.4 0.6 0.8 1.0 rc 10 20 30 40 50 1 1-rc -1 Fig. 1. 1/(1 −rc)−1 vs rc 0 0.2 0.4 0.6 0.8 1 0 500 1000 1500 2000 2500 rc N (c=50, ρ=0.8, Pim=0.2) µr=0.05 µr=0.1 µr=0.5 µr=1 µr=10 Fig. 2. rcvs Nfor c= 50, ρ= 0.8, Pim = 0.2 and µ= 1 From equation (4), the tail of the distribution e vjis geometrically distributed with parameter rc. In the stable state of the system described by CTMC e Y, the higher the value Nwe choose, the smaller the value of vNis. It is anticipated that the higher the chosen value of Nis the smaller the value of rcis if the system is in the stable state (see Figure 2). As one observes from Figure 1, the slope of the tangent line to the curve 1/(1 −rc)−1 increases as rcapproaches 1, and M(N)−Ntoo. Let rth be the selected upper limit of rc(the purpose of the selected upper limit is to search 19
for the initial value of Nto save the computational time), i.e., 0 < rc≤rth. Algorithm 3 To choose N 1: N←Nini 2: repeat 3: N←N+ 1 4: until lc+1(rth, N −1 + 1/(1 −rth)≤0 Algorithm 4 Proposed algorithm 1: Call algorithm 3 to choose N 2: step ←1 3: repeat 4: N←N+step 5: M0←N 6: k←0 7: repeat 8: k←k+ 1 9: Compute the root xcof lc+1(x, Mk−1) in the interval (0, rth] 10: Mk=N+rc 1−rc(using equation (12)) 11: iteration error =|Mk−Mk−1|/Mk−1 12: if lc+1(xc, Mk−1)> ǫrthen 13: k←0 14: N←N+ 1 15: M0←N 16: iteration error = 1 17: end if 18: until iteration error < ǫM 19: Compute ψcand R 20: Call Algorithm 2 21: Compute performance measures 22: Compute Converge 23: until Converge As stated before we also need to find the initial value of N. To save the computational time of the search for the initial value of N, we only check whether lc+1(x, M(N)) has a root in the interval (0, rth] instead of determining rc. Because lc+1(0, M(N)) = λholds, we have to examine whether lc+1(rth, M(N)) ≤ 0 in the first stage of our proposed solution. However, M(N) is not known in advance. To resolve this problem, Theorem 2 is applied to choose the initial value of N(see Algorithm 3), where lc+1(rth, N −1+1/(1−rth)) ≤0 is verified. Theorem 2 If rcis the eigenvalue of Q(x, M(N)) in the interval (0,1), and lc+1(rth, N −1 + 1/(1 −rth)) ≤0for rth >0, then rcis bounded by rth (i.e., 0< rc≤rth). 20
Proof. Assume that rth < rc, which follows lc+1(rth, M(N)) >0 because rcis the eigenvalue of Q(x, M(N)) in the interval (0,1) (i.e., lc+1(rc, M(N)) = 0) and lc+1(0, M(N)) = λ. Note that lc+1(x, y) is a monotonically decreasing function with respect to y for any constant x, 0 < x < 1, because •lc+1(0, y) = λholds, •lc+1(x, y) has a single root with respect to xin the interval (0,1) for any constant y, and •the higher the value of yis the smaller the root with respect to xis. We have N−1+1/(1−rth)< M(N) = N−1+1/(1−rc) from the assumption. Therefore, lc+1(rth, N −1 + 1/(1 −rth)) > lc+1(rth, M(N)). Using lc+1(rth, M(N)) >0, we obtain lc+1(rth, N −1 + 1/(1 −rth)) >0, which contradicts the given condition lc+1(rth, N −1 + 1/(1 −rth)) ≤0. Therefore, rth < rcdoes not hold, which yields that rcis bounded by rth.✷ Theorem 2 is used in our proposed Algorithm 4 to find the initial value of N. 3.5 The Convergence Criterion of a Proposed Algorithm We present the details of a proposed computational procedure in Algorithm 4 that integrates key results in Section 3.1-3.4. In Algorithm 4 two loops are applied after the initial choice of N. The inner loop is needed to find M(N) of a certain value N, while the outer loop is to tune Nto obtain the required accuracy of the estimation of performance measures. A notable feature of the proposed procedure compared to the HM2 algorithm is that the computation of M(N) does not require the computation of the steady-state probabilities. Furthermore, the algorithm is enhanced with the computation of threshold N. To define the convergence criterion Converge we made the following investigation: we compute the performance measures by running the core (between lines 4 and 19 of Algorithm 3) of our algorithm with fixing Nfor parameters µ= 1/180, µr= 0.01, Pim = 0.2, ǫM= 10−3and ǫr= 10−10 (Figures 3 and 4). From Figures 3 and 4 the performance measures converges as Ngrows. Furthermore, a high oscillation is observed in Figure 4 as well (see the e-companion for more numerical results with other settings of parameters). Therefore, the convergence criterion (to determine when the performance measures reach the “stable state”) is defined as the relative error between two moving averages 21
concerning a specific performance measure Converge := |PK i=Lχi/(K−L+ 1) −PK i=1 χi/K| PK i=Lχi/(K−L+ 1) < ǫp, where χi’s, i= 1,...,K, are the latest values of a chosen performance measure (e.g., Nret) determined so far by the algorithm, ǫpis the specified accuracy, K and Lare parameters. It is obvious that the minimum choice of Kand Lis K= 3 and L= 2. Furthermore, the higher the values of Kand Lare, the more time is needed, but the better the guarantee of convergence is ensured. In Table 1 we summarize the computational time of the algorithm for ρ= λ/(µc), µ= 1/180, µr= 0.01, Pim = 0.2. As observed the algorithm successfully stops in the stable region of performance measures (see Figures 3 and 4). From numerical results (see the e-companion for results with other settings of Kand L), the minimum choice of K= 3 and L= 2 can guarantee that the proposed algorithm steps over the “oscillation period”. Results in Table 1 and the e-companion empirically show that Pns can be set as the main performance measure in the convergence criterion Converge for determining all performance measures. 4 Computational Times of the Proposed Algorithm We plot the computational time versus cand Nin Figures 5 and 6, on a machine with IntelrCoreTM2 Duo T9400 2.53 GHz processor (note that the algorithm is implemented in Mathematica) for parameters ρ=λ/(µc) = 0.8, µ= 1, Pim = 0.2, ǫr= 10−10,ǫM= 10−3. In the curves the computational time of the HM2 algorithm [20] and the core (between lines 4 and 19 of Algorithm 3) of our algorithm. As observed the computational time of the original algorithm HM2 is a rapidly increasing function of c, while the computational complexity of the core of our algorithm is only O(c). In Figure 7, we plot the computational time of our algorithm versus ρand c. It is observed from Table 1 and Figure 7 that ρimpacts the convergence of the algorithms. The explanation behind this phenomenon is that the higher ρis the more likely the oscillation is (see Figure 4 and illustrations in the e-companion). Therefore, more computational time is needed to step over the oscillation. In other words, the algorithm needs more time to reach the convergence due to the oscillation phenomenon at large values of ρ. However, the 22
9.27e-008 9.271e-008 9.272e-008 9.273e-008 9.274e-008 9.275e-008 9.276e-008 9.277e-008 9.278e-008 9.279e-008 9.28e-008 20 40 60 80 100 Nret N 7.9e-009 7.92e-009 7.94e-009 7.96e-009 7.98e-009 8e-009 20 40 60 80 100 Pb N 7.84e-009 7.842e-009 7.844e-009 7.846e-009 7.848e-009 7.85e-009 7.852e-009 7.854e-009 7.856e-009 20 40 60 80 100 Pds N 9.6e-011 9.65e-011 9.7e-011 9.75e-011 9.8e-011 9.85e-011 9.9e-011 9.95e-011 1e-010 20 40 60 80 100 Pns N Fig. 3. Performance measures vs Nfor c= 50, ρ= 0.4, µ= 1/180, µr= 0.01, Pim = 0.2 23
74.4 74.405 74.41 74.415 74.42 74.425 74.43 74.435 74.44 74.445 74.45 80 100 120 140 160 180 Nret N 0.7521 0.75215 0.7522 0.75225 0.7523 0.75235 0.7524 0.75245 0.7525 80 100 120 140 160 180 Pb N 0.4618 0.46185 0.4619 0.46195 0.462 0.46205 0.4621 80 100 120 140 160 180 Pds N 0.2902 0.29025 0.2903 0.29035 0.2904 0.29045 0.2905 80 100 120 140 160 180 Pns N Fig. 4. Performance measures vs Nfor c= 50, ρ= 1.4, µ= 1/180, µr= 0.01, Pim = 0.2 24
Table 1 Nand computational time of the algorithm for K= 3, L = 2, ǫp= 10−5,ǫM= 10−3, ǫr= 10−10,ρ=λ/(µc), µ= 1/180, µr= 0.01, Pim = 0.2, rth = 0.95 c= 50 c= 100 c= 200 c= 500 c= 1000 ρ N Time (s) NTime (s) NTime (s) NTime (s) NTime (s) 0.4 Nret 17 0.889 6 0.765 8 1.435 39 16.318 13 13.121 Pb16 0.858 4 0.468 8 1.42 39 16.255 13 13.275 Pds 16 0.873 5 0.609 8 1.42 39 16.38 13 13.166 Pns 23 1.233 23 2.574 12 2.449 39 16.396 13 13.12 0.8 Nret 38 2.511 29 3.697 46 9.703 35 15.646 13 15.287 Pb38 2.574 32 4.57 41 8.408 35 15.647 12 13.744 Pds 38 2.511 32 4.446 41 8.377 35 15.709 11 11.762 Pns 47 3.182 32 4.617 46 9.751 36 17.035 40 52.355 1.0 Nret 44 2.855 76 7.504 117 14.133 206 35.491 363 82.619 Pb44 2.792 76 7.769 117 14.227 206 35.756 363 82.509 Pds 44 2.698 76 7.737 124 16.114 206 35.506 363 82.914 Pns 44 2.777 79 8.129 124 16.255 206 36.099 363 82.477 1.4 Nret 102 3.042 198 7.957 397 26.411 858 81.417 1679 264.001 Pb100 2.621 198 8.064 397 26.427 858 81.589 1679 266.122 Pds 106 3.463 198 8.003 397 26.318 858 81.339 1679 264.781 Pns 106 3.464 202 9.375 397 26.458 858 81.354 1679 266.013 computational time complexity is still of O(c) for a specific ρ. Remark. It is worth emphasizing that approaches belonging to the category “Approximations” (Domenech-Benlloch et al. [20]) produced unacceptable errors in most cases. Therefore, we do not focus on the comparison with these algorithms in this paper. The HM2 algorithm overcomes other approaches in the term of the accuracy. Our algorithm has the same accuracy as the HM2. However, it is much faster than the HM2 algorithm. We have shown that the computational time complexity of our algorithm is of O(c). To our best knowledge, we do not know that there is any other algorithm which has the computational time complexity of O(c) for the M/M/c retrial queue with impatient customers and has the same accuracy as of the HM2 algorithm (see Domenech-Benlloch et al. [20] and Do [3] for the overview of the latest algorithms). 25
0.87 0.872 0.874 0.876 0.878 0.88 20 40 60 80 100 Nret N 0.0331 0.0332 0.0333 0.0334 0.0335 0.0336 20 40 60 80 100 Pb N 0.0312 0.0314 0.0316 0.0318 0.032 0.0322 20 40 60 80 100 Pds N 0.0014 0.00145 0.0015 0.00155 0.0016 20 40 60 80 100 Pns N Fig. 2. Performance measures vs Nfor c= 50, ρ= 0.8, µ= 1/180, µr= 0.01, Pim = 0.2 3
14.32 14.33 14.34 14.35 14.36 14.37 14.38 14.39 14.4 14.41 14.42 20 40 60 80 100 120 Nret N 0.3324 0.3325 0.3326 0.3327 0.3328 0.3329 0.333 0.3331 0.3332 20 40 60 80 100 120 Pb N 0.286 0.2865 0.287 0.2875 0.288 20 40 60 80 100 120 Pds N 0.045 0.0455 0.046 0.0465 0.047 20 40 60 80 100 120 Pns N Fig. 3. Performance measures vs Nfor c= 50, ρ= 1.0, µ= 1/180, µr= 0.01, Pim = 0.2 4
74.4 74.405 74.41 74.415 74.42 74.425 74.43 74.435 74.44 74.445 74.45 80 100 120 140 160 180 Nret N 0.7521 0.75215 0.7522 0.75225 0.7523 0.75235 0.7524 0.75245 0.7525 80 100 120 140 160 180 Pb N 0.4618 0.46185 0.4619 0.46195 0.462 0.46205 0.4621 80 100 120 140 160 180 Pds N 0.2902 0.29025 0.2903 0.29035 0.2904 0.29045 0.2905 80 100 120 140 160 180 Pns N Fig. 4. Performance measures vs Nfor c= 50, ρ= 1.4, µ= 1/180, µr= 0.01, Pim = 0.2 5
4.16183e-069 4.16183e-069 4.16184e-069 4.16184e-069 4.16185e-069 4.16185e-069 4.16186e-069 4.16186e-069 4.16187e-069 20 40 60 80 100 Nret N 3.7275e-071 3.72755e-071 3.7276e-071 3.72765e-071 3.7277e-071 3.72775e-071 3.7278e-071 20 40 60 80 100 Pb N 3.723e-071 3.72305e-071 3.7231e-071 3.72315e-071 3.7232e-071 3.72325e-071 3.7233e-071 20 40 60 80 100 Pds N 4.48e-074 4.482e-074 4.484e-074 4.486e-074 4.488e-074 4.49e-074 4.492e-074 4.494e-074 4.496e-074 4.498e-074 4.5e-074 20 40 60 80 100 Pns N Fig. 5. Performance measures vs Nfor c= 500, ρ= 0.4, µ= 1/180, µr= 0.01, Pim = 0.2 6
4.11e-005 4.111e-005 4.112e-005 4.113e-005 4.114e-005 4.115e-005 4.116e-005 4.117e-005 4.118e-005 4.119e-005 4.12e-005 20 40 60 80 100 Nret N 1.81e-007 1.815e-007 1.82e-007 1.825e-007 1.83e-007 1.835e-007 1.84e-007 1.845e-007 1.85e-007 20 40 60 80 100 Pb N 1.81e-007 1.811e-007 1.812e-007 1.813e-007 1.814e-007 1.815e-007 1.816e-007 1.817e-007 1.818e-007 1.819e-007 1.82e-007 20 40 60 80 100 Pds N 7e-010 7.05e-010 7.1e-010 7.15e-010 7.2e-010 7.25e-010 7.3e-010 7.35e-010 7.4e-010 20 40 60 80 100 Pns N Fig. 6. Performance measures vs Nfor c= 500, ρ= 0.8, µ= 1/180, µr= 0.01, Pim = 0.2 7
65.1 65.2 65.3 65.4 65.5 65.6 65.7 100 150 200 250 Nret N 0.187 0.1875 0.188 0.1885 0.189 100 150 200 250 Pb N 0.1754 0.1756 0.1758 0.176 0.1762 0.1764 100 150 200 250 Pds N 0.0115 0.0116 0.0117 0.0118 0.0119 0.012 0.0121 0.0122 100 150 200 250 Pns N Fig. 7. Performance measures vs Nfor c= 500, ρ= 1.0, µ= 1/180, µr= 0.01, Pim = 0.2 8
737 737.5 738 738.5 739 740 760 780 800 820 840 860 880 Nret N 0.75325 0.753255 0.75326 0.753265 0.75327 0.753275 0.75328 0.753285 0.75329 740 760 780 800 820 840 860 880 Pb N 0.466 0.4665 0.467 0.4675 0.468 740 760 780 800 820 840 860 880 Pds N 0.285 0.2855 0.286 0.2865 0.287 740 760 780 800 820 840 860 880 Pns N Fig. 8. Performance measures vs Nfor c= 500, ρ= 1.4, µ= 1/180, µr= 0.01, Pim = 0.2 9
Table 1 Nand the computational time of the proposed algorithm for K= 3, L= 2 ǫp= 10−5,ǫM= 10−3,ǫr= 10−10,ρ=λ/(µc), µ= 1/180, µr= 0.01, Pim = 0.2, rth = 0.95 c= 50 c= 100 c= 200 c= 500 c= 1000 ρ N Time (s) NTime (s) NTime (s) NTime (s) NTime (s) 0.4 Nret 17 0.889 6 0.765 8 1.435 39 16.318 13 13.121 Pb16 0.858 4 0.468 8 1.42 39 16.255 13 13.275 Pds 16 0.873 5 0.609 8 1.42 39 16.38 13 13.166 Pns 23 1.233 23 2.574 12 2.449 39 16.396 13 13.12 0.8 Nret 38 2.511 29 3.697 46 9.703 35 15.646 13 15.287 Pb38 2.574 32 4.57 41 8.408 35 15.647 12 13.744 Pds 38 2.511 32 4.446 41 8.377 35 15.709 11 11.762 Pns 47 3.182 32 4.617 46 9.751 36 17.035 40 52.355 1.0 Nret 44 2.855 76 7.504 117 14.133 206 35.491 363 82.619 Pb44 2.792 76 7.769 117 14.227 206 35.756 363 82.509 Pds 44 2.698 76 7.737 124 16.114 206 35.506 363 82.914 Pns 44 2.777 79 8.129 124 16.255 206 36.099 363 82.477 1.4 Nret 102 3.042 198 7.957 397 26.411 858 81.417 1679 264.001 Pb100 2.621 198 8.064 397 26.427 858 81.589 1679 266.122 Pds 106 3.463 198 8.003 397 26.318 858 81.339 1679 264.781 Pns 106 3.464 202 9.375 397 26.458 858 81.354 1679 266.013 10
Table 2 Nand the computational time of the proposed algorithm for K= 4, L= 2, ǫp= 10−5,ǫM= 10−3,ǫr= 10−10,ρ=λ/(µc), µ= 1/180, µr= 0.01, Pim = 0.2, rth = 0.95 c= 50 c= 100 c= 200 c= 500 c= 1000 ρ N Time (s) NTime (s) NTime (s) NTime (s) NTime (s) 0.4 Nret 23 1.061 7 0.765 11 1.841 43 17.504 14 13.541 Pb17 0.78 5 0.484 11 1.809 43 17.487 14 13.588 Pds 17 0.749 5 0.499 11 1.779 43 17.254 14 13.462 Pns 23 1.06 36 3.213 14 2.559 43 17.316 14 13.447 0.8 Nret 47 2.808 30 3.713 53 10.639 36 15.694 15 16.614 Pb47 2.761 34 4.352 46 8.705 36 15.631 12 12.465 Pds 47 2.761 34 4.337 46 8.814 36 15.616 12 12.387 Pns 47 2.793 34 4.305 53 10.608 37 16.926 43 52.057 1.0 Nret 52 2.933 79 7.223 124 14.445 225 44.959 364 95.52 Pb52 2.932 79 7.332 124 14.617 225 44.85 364 95.02 Pds 52 2.933 79 7.363 124 14.633 225 44.975 364 94.584 Pns 61 3.698 83 7.971 138 18.362 225 44.788 364 94.833 1.4 Nret 106 3.042 202 8.409 404 28.735 863 98.421 1684 332.251 Pb102 2.62 202 8.362 404 28.657 863 98.515 1684 331.315 Pds 106 3.011 202 8.377 404 29.343 863 98.405 1684 330.16 Pns 106 3.042 202 8.362 404 28.735 863 98.078 1684 332.641 11
Table 3 Nand the computational time of the proposed algorithm for K= 5, L= 2, ǫp= 10−5,ǫM= 10−3,ǫr= 10−10,ρ=λ/(µc), µ= 1/180, µr= 0.01, Pim = 0.2, rth = 0.95 c= 50 c= 100 c= 200 c= 500 c= 1000 ρ N Time (s) NTime (s) NTime (s) NTime (s) NTime (s) 0.4 Nret 29 1.389 14 1.295 12 2.122 55 23.447 15 15.351 Pb23 1.045 6 0.655 12 2.122 55 23.447 15 15.21 Pds 23 1.061 6 0.624 12 2.138 55 23.416 15 15.288 Pns 29 1.341 37 3.447 18 3.245 55 23.588 15 15.288 0.8 Nret 49 2.932 32 3.993 56 11.715 37 17.316 16 18.642 Pb49 2.98 45 5.553 53 10.671 37 16.973 13 14.055 Pds 49 2.979 45 5.491 53 10.686 37 17.051 13 14.04 Pns 49 2.949 45 5.569 56 12.137 40 19.126 49 59.717 1.0 Nret 54 3.198 83 7.94 138 18.346 229 52.573 386 134.051 Pb54 3.198 83 8.003 138 18.58 229 52.65 386 132.102 Pds 54 3.104 83 8.018 138 18.626 229 52.276 386 133.583 Pns 54 3.182 83 8.05 140 20.03 229 52.261 386 133.162 1.4 Nret 111 3.573 204 9.547 407 32.947 872 123.319 1688 406.523 Pb106 3.01 204 9.454 407 32.76 872 123.287 1688 409.581 Pds 111 3.51 204 9.531 407 33.166 872 121.665 1688 411.032 Pns 111 3.573 204 9.563 407 32.838 872 122.321 1688 408.005 12