Double orbit finite retrial queues with priority customers and service interruptions

Applied Mathematics and Computation 253 (2015) 324–344

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Double orbit finite retrial queues with priority customers and service interruptions Madhu Jain a, Amita Bhagat a,⇑, Chandra Shekhar b a b

Department of Mathematics, Indian Institute of Technology Roorkee, Roorkee 247 667, India Department of Mathematics, Birla Institute of Technology and Science, Pilani, Rajasthan 333 031, India

a r t i c l e

i n f o

Keywords: Double orbits Unreliable server Finite capacity Threshold recovery policy Cost function Cellular radio network

a b s t r a c t The present study deals with the double orbit finite capacity retrial queues with unreliable server. The system facilitates the arrival of two types of customers known as priority and non priority customers and can hold a maximum of L priority customers and K non-priority customers as per its capacity. The priority customers are served prior to the non-priority customers. Moreover, the server is unreliable which may breakdown while servicing either priority or non-priority customer. The failed server is sent for repair following threshold recovery policy to become as good as earlier. Both transient as well as steady state analysis of the model has been done using by matrix method. Various performance measures including queue length, reliability metrics, long run probabilities, etc. have been obtained using various state probabilities. The application of the model to cellular radio network has been discussed. The cost function has been constructed and optimized using meta heuristic approach. The sensitivity analysis of various performance indices has been performed as an illustration. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction The formation of queues in front of railway reservation counters; doctor’s clinic, shopping malls etc. are day to day scenes which we encounter everyday in almost all the spheres of life throughout the world. Sometimes, due to congestion situations customers are deprived of service and hence decide to retry later on as soon as the congestion reduces or server becomes free. This state of reattempts for service gives rise to special queues known as retrial queues. These queues are the result of the fact that a customer who finds the server unavailable for service joins orbit which is a virtual pool of customers from where he tries again for the service as soon as he finds the server free. Enormous work has been done by a number of researchers in the direction of retrial queues and the related literature can be found in the survey articles by Yang and Templeton [1], Artalejo [2,3] Artalejo and Falin [4] and Artalejo [5]. The detailed account on retrial queues are given in the books on retrial queues by Falin and Templeton [6] and Artalejo and Coral [7]. In real life situations, the server seems to be unreliable due to which the service of the customers may be interrupted before its completion. It is worthwhile to have a look into the important works done on retrial queues. Li and Zhao [8] investigated queueing model with impatient customers in addition to constant retrial rate and unreliable server. An M/M/1 retrial queue with unreliable server has been investigated by Sherman and Kharoufeh [9] with infinite orbit capacity. Yang and Alfa

⇑ Corresponding author. E-mail addresses: [email protected] (M. Jain), [email protected] (A. Bhagat), [email protected] (C. Shekhar). http://dx.doi.org/10.1016/j.amc.2014.12.066 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

325

[10] studied a class of multi server queueing systems with server failures. Wang and Yang [11] used Quasi-Newton method to study the M/M/1/K queueing system with F-policy and unreliable server. Jain et al. [12] studied an unreliable M/M/2/K queueing system under (N,F) policy with multi optional phase repair. To assign priority to one type of jobs over other type of jobs is now a feature of daily routine activities. Now-a-days, while analyzing the queueing scenarios, almost every area seems to be affected by the priority service of one customer over another customer. Choi and Chang [13] analyzed various single server retrial queues with priority calls and many more features like geometric loss, feedback, etc. Artalejo et al. [14] made stationary analysis of a retrial queue with preemptive repeated attempts. Dimitriou [15] investigated a mixed priority retrial queue with negative arrivals, unreliable server and multiple vacations and explored stability conditions as well as the system state probabilities of the system. It usually happens in machining environment that the repairman is called upon only when a certain minimum number of faults are present in the system and are to be repaired as it is economical both in terms of time and expenses. Efrosinin and Semenova [16] studied M/M/1 system with unreliable server and threshold based repair policy. Efrosinin and Winkler [17] examined a Markovian retrial queueing system with non-reliable server. Jain and Bhagat [18] investigated a finite population retrial queueing system with threshold recovery and unreliable server with geometric arrivals and impatient customers. Purohit et al. [19] analyzed finite queue with threshold recovery and state dependent arrival rates. Jain and Bhagat [20] also studied a finite capacity retrial queueing system with finite capacity and impatient customers. The transient behavior of various queueing systems has also been studied by several researchers. Anisimov [21] studied averaging methods for the transient behavior in retrial queueing systems. Krinik et al. [22] examined the transient probabilities of single server retrial system using randomization solution. A new formula was developed for the transient solution of Erlang queueing system by Leonenko [23]. The present problem is framed in a very novel manner in the sense that the concept of double orbits is incorporated along with unreliable server and priority customers. No such work has been done as per authors’ information. The problem is equipped with retrial orbits having fixed capacity for both the priority and non priority customers. Moreover, the broken down server is repaired following threshold recovery for both priority as well as non priority customers. Non priority customers are served only if no priority customers are present in the system. The transient state solution of the model has been explored using matrix method. To highlight the practical importance of the model under consideration, numerical simulation has been done using an illustration. The cost function has been optimized to determine the optimal parameters. The rest of the paper is organized in the following manner. Section 2 deals with the description of the model with various assumptions and notations. The governing equations are established in Section 3. Various performance measures have been obtained in Section 4. An application on cellular mobile network has been developed in Section 5 and its numerical simulation has been done in Section 6. Section 7 deals with cost optimization of the system. Finally, the conclusions are drawn in Section 8.

2. Model description We consider an unreliable single server finite retrial queueing model with two types of customers; priority and non-priority customers with double orbits. The basic operation of the model can be described as: (i) Arrival and retrial process: Two class of customers namely priority and non-priority customers arrive in the system. The priority (non-priority) customers follow Poisson process with mean arrival rate k1 (k2). On finding the server idle or unavailable for the service, the customers on their arrival may either join the queue in the front of server or they can wait for their turn in their respective retrial orbits. The queue in the front of server can hold a maximum of L priority and K nonpriority customers. Further, we assume that the capacity of retrial orbit 1(orbit 2) is of L  1 (K  1) priority (non-priority) customers. The priority (non priority) customers retry for their service from their respective orbits; the retrial time is exponentially distributed with rate h1 (h2). The retrial process occurs only when single types of customers (either priority or non-priority) are present in the system. However, if both types of customers are present in the system then the priority customers are served like a classical queue and no retrial phenomenon takes place in such a case. (ii) Service process: All the customers are served following the discipline of first come first serve (FCFS); the priority and non priority customers are served according to exponential distribution with rates l1 and l2 ; respectively. The priority customers are served prior to the non-priority customers. The server after serving the last priority customer present in the queue automatically starts the servicing of non-priority customers waiting for the service. There is no retrial mechanism in this case. (iii) Breakdown state: The server is unreliable and may break down while serving the customers. The server failures occur in Poisson fashion with rates a1 (a2) while servicing priority (non priority) customers. The server failures occur while he is busy in serving either type of customer. No break down occurs during retrial and idle states. In case when the server breakdowns occurs, the service of the customer already in service is resumed and continued as soon as the repair process of the server is completed. Moreover, the system transition from repair states is allowed only into the busy states and not in the retrial states.

326

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

(iv) Repair process: The server may breakdown while serving either priority or non-priority customers. The broken down server is sent for repair so as to regain its functionality. The repair process is completed according to the threshold recovery policy: (a) If the server breakdown occurs while serving the priority customers, then the repair process starts if a sufficient number of priority customers say q1 (1 6 q1 6 L  1) are already present in the system. (b) If the server breakdown occurs while serving the non-priority customers, then the repair process starts if a threshold number of non-priority customers say q2 (1 6 q2 6 K  1) are available in the system. The various notations used to formulate the model are summarized below: k1 (k2): Arrival rate for priority (non priority) customers. l1, (l2): Service rate for priority (non priority) customers. h1 (h2): Retrial rate for priority (non priority) customers. a1 (a2): Breakdown rate for the server while servicing priority (non priority) customers. q1 (q2): Threshold limit on number of priority (non-priority) customers for the repair. b1 (b2): Repair rate for the down server broken during the service of priority (non priority) customers. 3. Mathematical formulation of model The retrial queueing model under consideration is markovian. In order to obtain the solution of the system we first develop mathematical model for the system using notations and assumptions discussed in previous Section 2. Chapman Kolmogorov equations are established to obtain the transient solution of the system. These are further explained in Sections 3.1 and 3.2. 3.1. Transient state probabilities Let ‘n’ and ‘m’ represent the number of priority and non priority customers present in the system (i.e. the customers in the orbit or queue including the customer with the server) at any time t, respectively. The transient state probabilities of the system are denoted as: P 0;0 ðtÞ = Probability of the server being in idle or inactive state with no customers in the system at time t. P n;m ðtÞ = Probability that the server is in busy state with n (0 6 n 6 L) priority and m (0 6 m 6 K) non-priority customers in the system at time t (except n = m = 0). Q ð1Þ n ðtÞ = Probability that there are n (1 6 n 6 L  1) priority customers in the retrial orbit 1 at time t. In this state, there is no service, however, the server is ready to work and the corresponding orbit is not empty. Q ð2Þ m ðtÞ = Probability that there are m (1 6 m 6 K  1) non priority customers in the retrial orbit 2 at time t, there is no service and the server is ready to work and the corresponding orbit is not empty. Rn;m ðtÞ = Probability that the server is under repair state with n (n P 0) priority and m (m P 0) non-priority customers in the system at time t (except n = m = 0). The steady state probabilities are given as: ð1Þ ð2Þ ð2Þ Pn;m ¼ limPn;m ðtÞ; Q ð1Þ n ðtÞ ¼ limQ n ðtÞ; Q m ðtÞ ¼ limQ m ðtÞ; Rn;m ¼ limRn;m ðtÞ t!1

Also, we denote

t!1

P 0n;m ðtÞ

¼

t!1

t!1

dP n;m ðtÞ . dt

3.2. Governing equations The Chapman–Kolmogorov equations for different states of the model (see Fig. 1) are constructed as written below. The following indicator functions are used for the formulation of the equations governing the model:

IA ¼



 IB ¼

1 if

ð0  n  q1  1Þ; ð0  m  K  1Þ

0 if

ðq1  n  L  1Þ; ð0  m  K  1Þ

1 if 0 if

ð0  m  q2  1Þ; ð0  n  L  1Þ ðq2  m  K  1Þ; ð0  n  L  1Þ

(i) Inactive state This state corresponds to the idle state of the server when neither priority nor non-priority customer is present in the system. The equation in this case is:

P00;0 ðtÞ ¼ ðk1 þ k2 ÞP0;0 ðtÞ þ l1 P1;0 ðtÞ þ l2 P0;1 ðtÞ

ð1Þ

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

327

Fig. 1. State transition diagram.

(ii) Busy state The present state corresponds to that state of the server when it is busy in rendering service to the customers. Depending on the number of customers present in the system, various cases can be framed as: (a) When there is no non-priority customer in the system ð1Þ

P01;0 ðtÞ ¼ ðk1 þ k2 þ a1 þ l1 ÞP 1;0 ðtÞ þ k1 P0;0 ðtÞ þ h1 Q 1 ðtÞ þ ð1  IA Þb1 R1;0 ðtÞ ð1Þ

P0n;0 ðtÞ ¼ ðk1 þ k2 þ a1 þ l1 ÞP n;0 ðtÞ þ k1 Pn1;0 ðtÞ þ h1 Q nð1Þ ðtÞ þ k1 Q n1 ðtÞ þ ð1  IA Þb1 Rn;0 ðtÞ; ð2  n  L  1Þ ð1Þ

P0L;0 ðtÞ ¼ ðk2 þ a1 þ l1 ÞPL;0 ðtÞ þ k1 PL1;0 ðtÞ þ k1 Q L1 ðtÞ þ b1 RL;0 ðtÞ

ð2Þ ð3Þ ð4Þ

(b) When both priority and non-priority customers are present in the system

P0n;m ðtÞ ¼ ðk1 þ k2 þ a1 þ l1 ÞPn;m ðtÞ þ k1 Pn1;m ðtÞ þ k2 Pn;m1 ðtÞ þ l1 P nþ1;m ðtÞ þ ð1  IA Þb1 Rn;m ðtÞ; ð1  n  L  1Þ; ð1  m  K  1Þ

ð5Þ

P0L;m ðtÞ ¼ ðk2 þ a1 þ l1 ÞPL;m ðtÞ þ k1 P L1;m ðtÞ þ k2 PL;m1 ðtÞ þ b1 RL;m ðtÞ; ð1  m  K  1Þ

ð6Þ

P0L;K ðtÞ ¼ ða1 þ l1 ÞPL;K ðtÞ þ k1 PL1;K ðtÞ þ k2 PL;K1 ðtÞ þ b1 RL;K ðtÞ

ð7Þ

328

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

P0n;K ðtÞ ¼ ðk1 þ a1 þ l1 ÞPn;K ðtÞ þ k1 Pn1;K ðtÞ þ k2 P n;K1 ðtÞ þ ð1  IA Þb1 RL;K ðtÞ þ l1 P nþ1;K ðtÞ; ð1  n  L  1Þ

ð8Þ

(c) When there is no priority customer in the system ð2Þ

P00;1 ðtÞ ¼ ðk1 þ k2 þ a2 þ l2 ÞP0;1 ðtÞ þ k2 P0;0 ðtÞ þ h2 Q 1 ðtÞ þ l1 P1;1 ðtÞ þ ð1  IB Þb2 R0;1 ðtÞ

ð9Þ

ð2Þ

P00;m ðtÞ ¼ ðk1 þ k2 þ a2 þ l2 ÞP0;m ðtÞ þ k2 P0;m1 ðtÞ þ k2 Q m1 ðtÞ þ h2 Q ð2Þ m ðtÞ þ l1 P 1;m ðtÞ þ ð1  I B Þb2 R0;m ðtÞ; ð2  m  K  1Þ

ð10Þ ð2Þ

P00;K ðtÞ ¼ ðk1 þ a2 þ l2 ÞP0;K ðtÞ þ k2 P0;K1 ðtÞ þ k2 Q K1 ðtÞ þ l1 P1;K ðtÞ þ b2 R0;K ðtÞ

ð11Þ

(iii) Repair state The server is unreliable and may break down while serving either the priority or non-priority customers. Therefore, present states correspond to the repair process of the broken down server in various situations as given below: (a) When the number of non priority customers = 0 and priority customers P 1

R01;0 ðtÞ ¼ ðk1 þ k2 ÞR1;0 ðtÞ þ a1 P1;0 ðtÞ

ð12Þ

R0n;0 ðtÞ ¼ ðk1 þ k2 þ ð1  IA Þb1 ÞRn;0 ðtÞ þ a1 Pn;0 ðtÞ þ k1 Rn1;0 ðtÞ; ð2  n  L  1Þ

ð13Þ

R0L;0 ðtÞ ¼ ðk2 þ b1 ÞRL;0 ðtÞ þ a1 PL;0 ðtÞ þ k1 RL1;0 ðtÞ

ð14Þ

(b) When the number of non priority customer P 1 and priority customers P 1

R0n;m ðtÞ ¼ ðk1 þ k2 þ ð1  IA Þb1 ÞRn;m ðtÞ þ a1 Pn;m ðtÞ þ k1 Rn1;m ðtÞ þ k2 Rn;m1 ðtÞ; ð1  n  L  1Þ; ð1  m  K  1Þ

ð15Þ

(c) When the number of non priority = K and priority customers P 1

R0n;K ðtÞ ¼ ðk1 þ ð1  IA Þb1 ÞRn;K ðtÞ þ a1 Pn;K ðtÞ þ k1 Rn1;K ðtÞ þ k2 Rn;K1 ðtÞ; ð1  n  L  1Þ

ð16Þ

(d) When the number of non priority customers P 1 and priority customers = L

R0L;m ðtÞ ¼ ðk2 þ ð1  IA Þb1 ÞRL;m ðtÞ þ a1 PL;m ðtÞ þ k1 RL1;m ðtÞ þ k2 RL;m1 ðtÞ; ð1  m  K  1Þ

ð17Þ

(e) When the number of priority customers = L and non priority customers = K

R0L;K ðtÞ ¼ b1 RL;K ðtÞ þ a1 P L;K ðtÞ þ k2 RL;K1 ðtÞ þ k1 RL1;K ðtÞ

ð18Þ

(f) When the number of priority customers = 0 and non priority customers P 1

R00;1 ðtÞ ¼ ðk1 þ k2 ÞR0;1 ðtÞ þ a2 P0;1 ðtÞ

ð19Þ

R00;m ðtÞ ¼ ðk1 þ k2 þ ð1  IA Þb2 ÞR0;m ðtÞ þ a2 P0;m ðtÞ þ k2 R0;m1 ðtÞ; ð2  m  K  1Þ

ð20Þ

(g) When the number of non priority customers = K and priority customers = 0

R00;K ðtÞ ¼ ðk1 þ b2 ÞR0;K ðtÞ þ a2 P0;K ðtÞ þ k2 R0;K1 ðtÞ

ð21Þ

(iv) Retrial orbit 1 (for priority customer) 0

ð1Þ Q ð1Þ for n ¼ 1; 2; 3; . . . ; L  1 n ðtÞ ¼ ðh1 þ k1 ÞQ n ðtÞ þ l1 P nþ1;0 ðtÞ

ð22Þ

(v) Retrial orbit 2 (for non priority customer) 0

ð2Þ Q ð2Þ for m ¼ 1; 2; 3; . . . ; K  1 m ðtÞ ¼ ðh2 þ k2 ÞQ m ðtÞ þ l2 P 0;mþ1 ðtÞ

ð23Þ

Normalizing condition is L1 K1 L X K L X K X X X X Q nð1Þ ðtÞ þ Q ð2Þ Pn;m ðtÞ þ Rn;m ðtÞ ¼ 1 m ðtÞ þ n¼1

m¼1

n¼0 m¼0

ð24Þ

n¼0 m¼0

4. Matrix method The present section deals with the transient solution of the system of equations by using the ‘Matrix method’. Matrix method is a well established technique to determine the solution of set of differential equations. At the initial stage, Laplace transforms of the equations are taken to convert them in differential free. The obtained new set of equations is then arranged

329

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

in the form of block matrix so as to obtain transient state probabilities in terms of eigenvalues of the determinant of coefficient matrix. The block matrix developed can be solved by a number of numerical techniques to obtain the eigen values of the coefficient matrix. The transient state probabilities are obtained in terms of the eigenvalues which can be further used to find the performance indices. Before proceeding further, we define transient state probabilities as:

PðtÞ ¼ ½P1 ðtÞ; P2 ðtÞ; P3 ðtÞ; . . . ; PKþ1 ðtÞ; PKþ2 ðtÞ; PKþ3 ðtÞT h

P1 ðtÞ ¼ P0;0 ðtÞ; Q 1ð2Þ ðtÞ; . . . ; Q ð2Þ K1 ðtÞ; P 0;1 ðtÞ; . . . ; P 0;K ðtÞ h

P2 ðtÞ ¼ R0;1 ðtÞ; . . . ; R0;K ðtÞ; Q 1ð1Þ ðtÞ; . . . ; Q ð1Þ L1 ðtÞ 

Plþ3 ðtÞ ¼ P 1;l ðtÞ; . . . ; P L;l ðtÞ; R1;l ðtÞ; . . . ; R1;l ðtÞ

iT ð2KÞð1Þ

iT ðKþL1Þð1Þ

T ð2LÞð1Þ

ð0  l  KÞ

For the sake of convenience, we define transient state probabilities in terms of

pðtÞ as-

P1 ðtÞ ¼ ½p1 ðtÞ; p2 ðtÞ; . . . ; p2K ðtÞTð2KÞð1Þ P2 ðtÞ ¼ ½p2Kþ1 ðtÞ; p2Kþ2 ðtÞ; . . . ; p3KþL1 ðtÞTðKþL1Þ1 

T

Plþ3 ðtÞ ¼ p3KþLþðl2LÞ ðtÞ; . . . ; p3KþL1þðlþ1Þ2LÞ ðtÞ

ð0  l  KÞ

ð2LÞð1Þ

ð25Þ

~ i ðsÞ. The Laplace transforms of pi ðtÞ is denoted by p ~ i ðsÞ; P ~ ðsÞ respectively. ~ i ðsÞ; P The Laplace transforms of pi ðtÞ; Pi ðtÞ; PðtÞ are denoted by p In terms of Laplace transform, the above transient state probabilities reduces to

h

~ ðsÞ ¼ P ~ 1 ðsÞ; P ~ 2 ðsÞ; P ~ 3 ðsÞ; . . . ; P ~ Kþ1 ðsÞ; P ~ Kþ2 ðsÞ; P ~ Kþ3 ðsÞ P

iT

ð26Þ

Noting that initially at time t = 0, the system is empty i.e.

P0;0 ð0Þ ¼ 1; Pn;m ð0Þ ¼ 0 8n – 0; m – 0

ð27Þ

The initial vector can be defined as,

Pð0Þ ¼ ½1; 0; 0; . . . ; 0; 0; 0; 0; . . . ; 0; 0; 0; 0; . . . ; 0ð2KLþ3Kþ3L1Þ1

ð28Þ

After taking their Laplace transforms, the set of differential equations (1)–(23) can be written in matrix form as:

~ ðsÞ ¼ Pð0Þ AðsÞP

ð29Þ

where,

2

A0

6 D 6 0 6 6 G 6 6 6 E1 6 6 AðsÞ ¼ 6 E2 6 6 E3 6 6 .. 6 . 6 6 4 EK1 EK  A0 ¼

B0 ¼



C

B

D

E

0 0 F 0

B0

B1

A1

G1

D1

A2

F1

D2

C1

C2

C3

A3 D2 ..

A3 ..

D2

.

FK1 FK

. D2



 A1 ¼

ð2KÞ

 ðKþL1Þð2KÞ

CK

A3

F2 F3 .. .

;

. . . CK1

H 0 0

; B1 ¼

I 



 ;

A2 ¼

ðKþL1Þ

J1 0

A3 D2

0 0

 ð2KÞð2LÞ

7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5

A4

M0

M

M1

U

; G1 ¼

3



ð2LKþ3K1Þð2LKþ3K1Þ



 ;

ð2LÞ

0 0 Z 0

A3 ¼

S

M

M1

U

 ðKþL1Þð2KÞ

; G¼



 A4 ¼

; ð2LÞ



T1 0

0 0

V

M

M1

W

 ð2LÞð2KÞ

 ð2LÞ

330

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

 D0 ¼  Cl ¼

0 R0 0

0

0 0 Jl

0



 ; D1 ¼

0

ðKþL1Þð2KÞ



 ; El ¼

ð2KÞð2LÞ

0 Tl 0

0 K1

0

0



 ; D2 ¼ ð2LÞðKþL1Þ



 ; Fl ¼ ð2LÞð2KÞ

0

0

Tl

0

R

0

0

R

 ð2LÞ

 l ¼ 1; 2; . . . K ð2LÞðKþL1Þ

The various sub matrices can be further defined as:

B ¼ diag:½l2 ðKKÞ ; R0 ¼ diag:½a2 ðKKÞ ; I ¼ diag:½ðs þ h1 þ k1 ÞðL1ÞðL1Þ ; M1 ¼ diag:½a1 ðLLÞ ; R ¼ diag:½k2 ðLLÞ Tl ¼ ½t ij ðLKÞ  k1 ; for i ¼ 1; j ¼ l tij ¼ 0; otherwise K1 ¼ ½kij ðLðL1ÞÞ  h1 ; for i ¼ j kij ¼ k1 ; for i ¼ j þ 1 and 1  j  L  1 M ¼ ½mij ðLLÞ  b1 ; for i ¼ j and q1  j  L mij ¼ 0; otherwise

F ¼ ½f ij ðKKÞ  b2 ; for i ¼ j and q2  j  K f ij ¼ 0; otherwise Jl ¼ ½jij ðKLÞ  l1 ; for i ¼ l and j ¼ 1 jij ¼ 0; otherwise

C ¼ ½aij ðKKÞ  ðs þ k1 þ k2 Þ; for i ¼ j ¼ 1 aij ¼ ðs þ h2 þ k2 Þ; for i ¼ j and 2  i  K

D ¼ ½dij ðKKÞ  k2 ; for i ¼ j; 1  i  K dij ¼ h2 ; for j ¼ i þ 1 and 1  i  K  1

Z ¼ ½zij ðL1ÞL  l1 ; for 1  i  L  1 and j ¼ i þ 1 zij ¼ 0; otherwise

E ¼ ½eij ðKKÞ 8 > < ðs þ k1 þ k2 þ a2 þ l2 Þ; for i ¼ j and 1  i  K  1 eij ¼ ðs þ k2 þ a2 þ l2 Þ; for i ¼ j and i ¼ K > : k2 ; for i ¼ j þ 1 and 1  j  K  1 M0 ¼ ½@ ij ðLLÞ 8 > < ðs þ k1 þ k2 þ a1 þ l1 Þ; for i ¼ j and 1  i  L  1 @ ij ¼ ðs þ k2 þ a1 þ l1 Þ; for i ¼ j and i ¼ L > : k1 ; for i ¼ j þ 1 and 1  j  L  1 U ¼ ½uij ðLLÞ 8 ðs þ k1 þ k2 Þ; for i ¼ j and 1  i  q1  1 > > > < ðs þ k1 þ k2 þ b1 Þ; for i ¼ j and q1  i  L  1 uij ¼ > ðs þ k2 þ b1 Þ; for i ¼ j ¼ L > > : k1 ; for i ¼ j þ 1 and 1  j  L  1 S ¼ ½sij ðLLÞ 8 ðs þ k1 þ k2 þ a1 þ l1 Þ; for i ¼ j and 1  i  L  1 > > > < ðs þ k þ a þ l Þ; for i ¼ j and i ¼ L 2 1 1 sij ¼ > k ; for i ¼ j þ 1 and 1jL1 1 > > : l1 ; for j ¼ i þ 1 and 1  i  L  1 H ¼ ½hij ðKKÞ 8 ðs þ k1 þ k2 Þ; for i ¼ j and 1  i  q2  1 > > > < ðs þ k þ k þ b Þ; for i ¼ j and q  i  K  1 1 2 2 2 hij ¼ > ðs þ k þ k þ b Þ; for i ¼ j ¼ K 1 2 > 2 > : k2 ; for i ¼ j þ 1 and 1  j  K  1

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

331

W ¼ ½wij ðLLÞ 8 ðs þ k1 Þ; for i ¼ j and 1  i  q1  1 > > > < ðs þ k1 þ b1 Þ; for i ¼ j and q1  i  L  1 wij ¼ > ðs þ b1 Þ; for i ¼ j ¼ L > > : k1 ; for i ¼ j þ 1 and 1  j  L  1 V ¼ ½v ij ðLLÞ 8 ðs þ k1 þ a1 þ l1 Þ; for i ¼ j and 1  i  L  1 > > > < ðs þ a1 þ l1 Þ; for i ¼ j ¼ L v ij ¼ > k1 ; for i ¼ j þ 1 and 1  j  L  1 > > : l1 ; for j ¼ i þ 1 and 1  i  L  1 Now, we aim to find out the probability of the system at any time t. At the initial stage, we use Cramer’s rule so as to determine the transient probabilities of the server at different states on matrix AðsÞ: From Eq. (29), we obtain (cf. Jain et al. [12]):

p~ i ðsÞ ¼

det½Ai ðsÞ ; ði ¼ 1; 2; 3; 4; . . . ; ð2LK þ 3L þ 3K  1ÞÞ det½AðsÞ

ð30Þ

where, det½AðsÞ is the determinant of the matrix A(s) and det½Ai ðsÞ is the determinant of matrix which has been obtained by replacing the respective ith column vector, ði ¼ 1; 2; 3; 4; . . . ; ð2LK þ 3L þ 3K  1ÞÞ of AðsÞ with initial vector Pð0Þ. Now, in order to obtain the explicit expression for Eq. (30), we proceed as follows: It is clear that s = 0 is a root of det½AðsÞ = 0. Now substituting s = d, we obtain

AðdÞ ¼ ðB  dIÞ

ð31Þ

where, B = Að0Þ is a square matrix of order ð2LK þ 3L þ 3K  1Þ and I is the identity matrix of order ð2LK þ 3L þ 3K  1Þ: Now using Eqs. (29) and (31), we obtain

~ ðsÞ ¼ ðB  dIÞP ~ ðsÞ ¼ Pð0Þ AðdÞP

ð32Þ

Now, we want to find out other distinct eigenvalues dx ðdx – 0; x ¼ 1; 2; . . . ; ð2LK þ 3L þ 3K  2ÞÞ of the matrix B  dI: For this purpose we equate its determinant equals to zero. The eigenvalues are either, real or complex. We assume that there are x real and y pairs of distinct conjugate complex eigenvalues which we denote by

d1 ; d2 ; . . . ; dx and ðdxþ1 ; dxþ1 Þ; ðdxþ2 ; dxþ2 Þ; . . . ; ðdxþy ; dxþy Þ; respectively Moreover, x þ 2y ¼ ð2LK þ 3L þ 3K  2Þ, where x is the number of real eigenvalues (excluding zero) and y represents the number of pairs of conjugate complex eigenvalues. Thus, we have

      det½AðsÞ ¼ s Pxk¼1 ðs þ dk Þ Pyk¼1 ðs þ dxþk Þðs þ dxþk Þ ¼ s Pxk¼1 ðs þ dk Þ Pyk¼1 fs2 þ ðdxþk þ dxþk Þs þ dxþk dxþk g

ð33Þ

Using Eqs. (30) and (33), we get

det½Ai ðsÞ   ; ði ¼ 1; 2; 3; 4; . . . ; ð2LK þ 3L þ 3K  1ÞÞ s Pxk¼1 ðs þ dk Þ Pyk¼1 fs2 þ ðdxþk þ dxþk Þs þ dxþk dxþk g

p~ i ðsÞ ¼ 

ð34Þ

On expanding by partial fractions, we get

p~ i ðsÞ ¼

y x X a0 X am bm s þ cm þ þ ; ði ¼ 1; 2; 3; 4; . . . ; ð2LK þ 3L þ 3K  1ÞÞ s m¼1 s þ dm m¼1 s2 þ ðdxþm þ dxþm Þs þ dxþm dxþm

ð35Þ

where a0 and am ðm ¼ 1; 2; . . . ; xÞ are real numbers and are obtained as

a0 ¼ 

am ¼

det½Ai ð0Þ 

Pxk¼1 dk Pyk¼1 dxþk dxþk



x k¼1k–m ðdk

ðdm Þ P



det½Ai ðdm Þ oi ; m ¼ 1; 2; . . . ; x h y n  dm Þ Pk¼1 ðdm Þ2 þ ðdxþk þ dxþk Þðdm Þ þ dxþk dxþk

ð36Þ

ð37Þ

332

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

bm ðdxþm Þ þ cm ¼

det½Ai ðdxþm Þ n oi ; m  h ðdxþm Þ Pxk¼1 ðdk  dxþm Þ Pyk¼1k–m ðdxþm Þ2 þ ðdxþk þ dxþk Þðdxþm Þ þ dxþk dxþk

¼ 1; 2; . . . ; y

ð38Þ

On taking inverse Laplace transformation of Eq. (35), the probability of the system state at any time t, are given by

pi ðtÞ ¼ a0 þ

 y  x X X cm  bm v m v m t am edm t þ bm ev m t cosðwm tÞ þ e sinðwm tÞ ; ð1  i  ð2LK þ 3L þ 3K  1ÞÞ wm m¼1 m¼1

where, a0 ; am ; bm ; cm ;

v m;

ð39Þ

wm are real numbers.

5. Matrix recursive approach In the previous section we used Matrix method to obtain transient state probabilities of the system in terms of eigenvalues of the determinants. Now, to obtain the steady state probabilities of the system we use ‘Matrix Recursive Approach’. This method is widely used to deal with various queueing models that demands for exact steady state solutions. The steady state probabilities are given as:

Pn;m ¼ limPn;m ðtÞ;

Q nð1Þ ðtÞ ¼ limQ nð1Þ ðtÞ;

Rn;m ¼ limRn;m ðtÞ;

Pn ¼ limPn ðtÞ

t!1

t!1

t!1

ð2Þ ð2Þ Qm ðtÞ ¼ limQ m ðtÞ t!1

ð40Þ

t!1

Eq. (29) can be written as:

2

A0

B0

6 D 6 0 6 6 G 6 6 6 E1 6 6 6 E2 6 6 E 6 3 6 6 .. 6 . 6 6 4 EK1 EK

B1

A1

G1

D1

A2

F1

D2

C1

C2

C3

A3

3 2 P1 ð0Þ 3 CK 2 3 P 2 7 7 6 76 7 6 P2 ð0Þ 7 76 P3 7 6 7 76 7 6 P3 ð0Þ 7 76 P 7 6 7 76 7 6 4 7 76 7 6 P4 ð0Þ 7 76 P5 7 6 7 76 7 6 P ð0Þ 7 76 7 7¼6 5 76 .. 7 6 7 76 . 7 6 7 . 76 7 7 6 .. 76 7 7 6 76 PKþ1 7 6 7 76 7 6 PKþ1 ð0Þ 7 76 P 7 7 6 74 Kþ2 5 6 7 5 4 PKþ2 ð0Þ 5

D2

A4

. . . CK1

A3

F2

D2

F3 .. .

A3 D2 ..

.

A3 ..

FK1

. D2

FK

PKþ3

PKþ3 ð0Þ

Therefore we have the following set of equations:

A0 P1 þ B0 P2 þ B1 P3 þ C1 P4 þ    þ CK1 PK1 þ CK PK ¼ 1

ð41Þ

D0 P1 þ A1 P2 þ G1 P3 ¼ 0

ð42Þ

GP1 þ D1 P2 þ A2 P3 ¼ 0 E1 P1 þ F1 P2 þ D2 P3 þ A3 P4 ¼ 0

ð43Þ ð44Þ

E2 P1 þ F2 P2 þ D2 P4 þ A3 P5 ¼ 0

ð45Þ

E3 P1 þ F3 P2 þ D2 P5 þ A3 P6 ¼ 0 .. .

ð46Þ

EK1 P1 þ FK1 P2 þ D2 PKþ1 þ A3 PKþ2 ¼ 0

ð47Þ

EK P1 þ FK P2 þ D2 PKþ2 þ A4 PKþ3 ¼ 0

ð48Þ

On solving (42) and (43), we get 1 P3 ¼ G1 1 fD0 P1 þ A1 P2 g and P3 ¼ A2 fGP1 þ D1 P2 g

n

1 P2 ¼ G1 1 A1  A2 D1

o1 n

o 1 A1 2 G  G1 D0 P1

Using the value of P2 in P3 ¼ G1 1 fD0 P1 þ A1 P2 g; we get



n

1 P3 ¼ G1 D0 þ G1 1 1 A1  A2 D1

o1 n

1 A1 2 G  G1 D0

o

P1

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

333

Proceeding in a similar pattern, we can obtain a general result for Eqs. (44)–(47) in the following form:

Pi ¼ A1 3

( i3

X

Ei2r ðD2 A1 3 Þ

r1



P1 þ

r¼1

i3 X

)

r1 i3 1 i4 Fi2r D2 A1 P þ ðD Þ ðA Þ P 2 2 3 ; 4  i  K þ2 3 3

r¼1

A1 4 fEK

Eq. (48) yields, PKþ3 ¼ P1 þ F K P2 þ D2 PKþ2 g. By using the values of unknowns from P2 to PKþ2 in Eq. (41), i.e.

A0 P1 þ B0 P2 þ B1 P3 þ

K X Ci Piþ3 ¼ 1; we can get P1 i¼1

6. Performance measures The validity of any retrial queueing model can be best deciphered in terms of its performance indices. Various indices namely average queue length, system state probabilities, throughput etc. can be determined so as to judge the efficiency of the system. Some of the performance measures are as follows: 6.1. Server state probabilities The probabilities of the different states of the server are important in deciding the efficiency and other metrics of the system. The various server state probabilities for the different states of the server at time t are established as follows: (a) Busy state: The probability that the server being busy in providing service to the customers at time t, is

PB ðtÞ ¼

L X K L K X X X P n;m ðtÞ þ Pn;0 ðtÞ þ P0;m ðtÞ n¼1 m¼1

n¼1

ð49Þ

m¼1

(b) Broken down state: The probability that the server being in broken down state at time t, is

Pf ðtÞ ¼

L X K X

Rn;m ðtÞ

ð50Þ

n¼1 m¼1

(c) Repair state: (i) The probability that the server starts the repair of server failed while servicing non-priority customers at time t, is K X

PR1 ðtÞ ¼

R0;m ðtÞ

ð51Þ

m¼q2

(ii) The probability that the server starts the repair of the server failed while servicing priority customers at time t, is

PR2 ðtÞ ¼

L X K X Rn;m ðtÞ

ð52Þ

n¼q1 m¼0

(iii) The probability that the server starts the repair of failed server at time t, is

PR ðtÞ ¼ P R1 ðtÞ þ P R2 ðtÞ

ð53Þ

(iv) The probability that the server is in broken down state but repair is not started at time t

PR3 ðtÞ ¼

qX K 1 1X

qX 2 1

n¼1 m¼0

m¼1

Rn;m ðtÞ þ

R0;m ðtÞ

ð54Þ

(d) Retrial state: (i) The probability that the priority customer retry for the service at time t, is

PRE1 ðtÞ ¼

L1 X Q ð1Þ n ðtÞ n¼1

ð55Þ

334

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

(ii) The probability that the non-priority customer retry for the service at time t, is

PRE2 ðtÞ ¼

K1 X

ð2Þ Qm ðtÞ

ð56Þ

m¼1

6.2. Queueing indices The computation of queueing indices of any system is the most significant and promising measure to upgrade any system. It really helps the system designers for better management of delay situations and efficient functioning of the queueing systems. (a) Queue length: The assessment of queue length is the primary objective of any queueing model. Here we give some indices which are related to queue length of the concerned model. (i) Expected number of priority customers in the system at time t

EfN1 ðtÞg ¼

L1 L X K X X nQ ð1Þ nðPn;m ðtÞ þ Rn;m ðtÞÞ n ðtÞ þ n¼1

ð57Þ

n¼1 m¼0

(ii) Expected number of non-priority customers in the system at time t

EfN2 ðtÞg ¼

K 1 X

mQ ð2Þ m ðtÞ þ

L X K X mðPn;m ðtÞ þ Rn;m ðtÞÞ

m¼1

ð58Þ

n¼0 m¼1

(iii) Expected number of customers in the system at time t

EfNðtÞg ¼ EfN1 ðtÞg þ EfN2 ðtÞg

ð59Þ

(iv) Expected number of customers waiting in the retrial orbits at any time t

EfNr ðtÞg ¼

L1 K1 X X ð2Þ nQ ð1Þ mQ m ðtÞ n ðtÞ þ n¼1

ð60Þ

m¼1

(v) Expected number of customers in the breakdown state at any time t

EfNf ðtÞg ¼

L X K L X K X X nRn;m ðtÞ þ mRn;m ðtÞ n¼1 m¼0

ð61Þ

n¼0 m¼1

(b) Throughput: Throughput can be considered as the average rate of successful services rendered to the customers by the server in a queueing system. It can be mathematically expressed as:

TPðtÞ ¼ l1

L X K K X X Pn;m ðtÞ þ l2 P 0;m ðtÞ n¼1 m¼0

ð62Þ

m¼1

(c) Carried load: The carried load at time t is:

" # " # L1 L X K L X K K1 L X K L X K X X X X X X ð1Þ ð2Þ C L ðtÞ ¼ k1 Q n ðtÞ þ Pn;m ðtÞ þ Rn;m ðtÞ þ k2 Q m ðtÞ þ Pn;m ðtÞ þ Rn;m ðtÞ n¼1

n¼1 m¼0

n¼1 m¼1

m¼1

n¼0 m¼1

ð63Þ

n¼1 m¼1

6.3. Reliability measures Like queueing indices and server state probabilities, reliability indices also play a major role in increasing the availability and efficiency of the concerned unreliable server queueing system. These indices can be further used to improve the system during design and development phases at time t. (i) Availability of the server at any time t It measures the probability of the server being available in the system and is given by

Av ðtÞ ¼

L1 K1 L X K X X X Q nð1Þ ðtÞ þ Q ð2Þ Pn;m m ðtÞ þ n¼1

m¼1

n¼0 m¼0

ð64Þ

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

335

(ii) Failure frequency The rate of failure of the server at different states is used to obtain the failure frequency of the system. Thus, the failure frequency at any time t is expressed as:

F f ðtÞ ¼ a1

L X K K X X Pn;m ðtÞ þ a2 P0;m ðtÞ n¼1 m¼0

ð65Þ

m¼1

7. Application of the model in cellular radio network We illustrate the real life application of queueing model under consideration in cellular radio network wherein the whole geographical area is divided in cells and in each cell there is a base station. We consider the radio transmission via base station by considering a single channel in the microcell to serve the incoming calls which are generated in Poisson fashion. The traffic flow in cellular radio network is clearly represented by means of Fig. 2. Further, we assume that before arrival of any call in the system, ‘n’ handoff calls and ‘m’ new calls are already present in the system. The arriving calls are of two types i.e. (i) new calls which are assumed to be originated in the cell in the coverage area of the cell and (ii) handover calls which are transferred from the neighboring cell due to mobility of the subscribers. A priority is given to the handoff calls as once connection is established and the call should be continued till completion by reserving some channels called guard channels for them but not at the cost of new calls. New call requests are allowed to be queued in the buffer whenever no channel is available at the arrival instant; in such a case the handover calls are treated as the blocked call. The interruption in the transmission due to unavailability of channel is managed by providing the buffer of capacity L (K) for handover (new) attempts. If any handover call is present in the cell, the new call is not served. If the channel is busy, the new as well as handover attempts wait in the orbits; the calls can retry for their attempts from the orbits. The channel is subject to breakdown and repair. The repairing of the channel is permissible when a pre-specified number of calls are accumulated in the buffer of coverage area. Here, we consider that handover calls arrive with rate k1 and new calls arrive with arrival rate k2 : The handover (new) calls are served with rate l1 (l2 ). Excluding the call being served at the moment by the channel, extra handover and new calls present in the queue retry for the service with retrial rates h1 (h2 ). A handoff call can retry for the service only if no new calls are present in the system, i.e. retrials cannot be made by handoff calls at the loss or cost of new (non-priority) calls. Hence, we can say that cellular radio network can be considered as the direct implication of our retrial model. 8. Numerical simulation In this section, we provide the numerical simulation by taking the illustration of the cellular mobile network as explained in previous Section 5. It is worth-mentioning that the computation of transient probabilities by the suggested matrix method

Fig. 2. Flow chart of traffic in cellular radio network.

336

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

will portray more realistic scenario of concerned retrial model for which closed form analytical results are not possible to derive. The numerical experiment is performed by developing code in MATLAB software to validate the tractability and applicability of the proposed computational approach in the concerned retrial model. We divide this section in various categories so as to study the effect of various parameters on the sensitivity of the system. To illustrate the system numerically, the values for default parameters are considered as

k1 ¼ 5; k2 ¼ 6; h1 ¼ 0:5; h2 ¼ 1;

l1 ¼ 8; l2 ¼ 7; b1 ¼ 1; b2 ¼ 2; a1 ¼ 0:5; a2 ¼ 0:5

8.1. Queue length of the system The single channel (server) can serve only one call at a time, therefore rest of the calls either handoff or new form queue in the system. The queue length of either type of calls is affected by various parameters of the system namely service rate of the server to attend the calls, arrival rate of calls and so on. Illustration 1: We consider the set of default parameters and vary q1 from 1 to 5 at a time range of 1 to 10 units and find corresponding queue length of priority customers. To minimize the congestion in the system, we want to find out the optimal threshold recovery parameter q1 which minimizes the queue length of priority customers (handover calls). Table 1 below displays the value of queue length of priority customers (handover calls) i.e. E[N1(t)] corresponding to different threshold parameters q1 and time t. It is clear from the table that the minimum queue length is obtained for q1 = 2. Fig. 3 diagrammatically represents the various values of queue length obtained in Table 1. It is very obvious that the number of calls increases with an increase in time t. However, E[N1(t)] decreases from q1 = 1 to q1 = 2 and then increases up to q1 = 5. It is very clear from the data depicted in above table that q1 = 2 seem to be the optimal threshold recovery point where minimum queue length is observed. Hence, with this set of default parameters, a minimum queue length for handover calls can be obtained if recovery of failed server starts corresponding to q1 = 2. This may help in reducing the congestion in the system. Figs. 4a and 4b depict the behavior of queue length of handoff calls i.e. priority customers w.r.t. arrival rate ðk1 Þ and breakdown rate ða1 Þ on E[N1(t)]. An increase in arrival rate (breakdown rate) increases the accumulation of handoff calls in the system (make the system more prone to failures). It reveals that the system can run smoothly if breakdown rate can be lowered. Figs. 5a–5c exhibit the variation in the queue length of both handoff calls (priority customers) and new calls (non-priority customers) with time t and other parameters. The variation in the queue length of new incoming calls (non-priority customers) with time t and parameters q2 and K have been demonstrated in Figs 5b–5c by means of surface graphs. The queue length E[N2(t)] of the incoming calls (non-priority customers) increases with an increase in time t as well as with q2 and K.

Table 1 Variation in E[N1(t)] with q1 and t.

Fig. 3. Effect of q1 on E[N1(t)] with time t.

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

E[N1(t)]

4 3 λ1=4 λ1=4.5 λ1=5 λ1=5.5 λ1=6

2 1 0 0

2

4 t

6

8

10

Fig. 4a. Effect of k1 on E[N1(t)] with t.

5 E[N1(t)]

4 3 α1=0.5 α1=1 α1=1.5 α1=2

2 1 0 0

2

4

t 6

8

Fig. 4b. Effect of a1 on E[N1(t)] with t.

Fig. 5a. Effect of L on E[N1(t)] with time.

Fig. 5b. Effect of q2 on E[N2(t)] with time t.

Fig. 5c. Effect of K on E[N2(t)] with time t.

10

337

338

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

8.2. Throughput Throughput can be defined as the amount of effective service rendered per unit time. Figs, 6a and 6b display the effect of service rate (l1) and arrival rate (k1) respectively, on the throughput TP(t) at time t. TP(t) is maximum in both the graphs when l1 = k1. However, in Fig. 6a we can observe that at t = 2 units and onwards, TP(t) shows a steady state behavior even on increasing the service rate l1. In Fig. 6b, TP(t) increases with the growth of k1; this is due to the fact that growth in the number of new calls automatically increases TP(t). Concluding we say that TP(t) can be enhanced by increasing the service rates for either type of calls. Also, higher TP(t) means higher efficiency of the system. 8.3. Reliability indices Similar to queueing indices, reliability measures facilitates the valuable information to justify and validate the effectiveness of queueing model of cellular radio network. Tables 2 and 3 present the effect of various parameters on the reliability indices; namely availability and failure frequency on the system with time t. Both availability Av(t) and failure frequency ðF f ðtÞÞ of the system decrease as time t grows. From data depicted in Table 2, it is clear that Av(t) increases but F f ðtÞF decreases with the growth of service rate l1 from 6 units to 8 units. This is due to the fact that an increase in service rate

7 6

TP(t)

5 4 3

μ1=5 μ1=6 μ1=7 μ1=8

2 1 0

0

1

2

3

t

4

5

Fig. 6a. Effect of l1 on TP(t) with time t.

6 5

TP(t)

4 3 λ1=5 λ1=6 λ1=7 λ1=8

2 1 0

0

1

2

t

3

4

5

Fig. 6b. Effect of k1 on TP(t) with time t.

Table 2 Effect of k1 on Av(t) and F f ðtÞ. k1

t

Av ðtÞ

F f ðtÞ

l1 = 6

l1 = 7

l1 = 8

l1 = 6

l1 = 7

l1 = 8

3

1 3 5

0.7837 0.6894 0.6852

0.7901 0.6974 0.6936

0.7966 0.7041 0.7009

0.2718 0.2490 0.2545

0.2680 0.2404 0.2445

0.2628 0.2338 0.2369

3.5

1 3 5

0.7802 0.6884 0.6832

0.7860 0.6962 0.6915

0.7915 0.7031 0.6990

0.2796 0.2613 0.2663

0.2780 0.2514 0.2550

0.2723 0.2436 0.2462

4

1 3 5

0.7743 0.6869 0.6810

0.7830 0.6946 0.6891

0.1059 0.6872 0.6963

0.2762 0.2727 0.2771

0.2858 0.2618 0.2649

0.2219 0.2314 0.2547

339

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344 Table 3 Effect of k2 on Av(t) and F f ðtÞ. t

k2

Av(t)

F f ðtÞ

b1 = 0.8

b1 = 1.2

b1 = 1.4

b1 = 0.8

b1 = 1.2

b1 = 1.4

4

1 3 5

0.7913 0.6781 0.6685

0.8163 0.7369 0.7361

0.8415 0.7575 0.7585

0.1194 0.2427 0.2431

0.2766 0.2666 0.2674

0.3044 0.2744 0.2751

5

1 3 5

0.7869 0.6710 0.6593

0.7987 0.7299 0.7278

0.8046 0.7509 0.7505

0.2790 0.2505 0.2505

0.2762 0.2746 0.2763

0.2770 0.2829 0.2845

6

1 3 5

0.7765 0.6640 0.6516

0.7909 0.7236 0.7210

0.7968 0.7449 0.7440

0.2951 0.2564 0.2563

0.3059 0.2813 0.2833

0.3100 0.2900 0.2920

of the server makes the server more available for the service. Further, Av(t) ðF f ðtÞÞ exhibit decreasing (increasing) behavior with the growth of arrival rates either k1 or k2. An increase in arrival rate makes the server busier with the customers and thus it becomes less available for new incoming calls. Based on numerical results for reliability indices, we infer that the availability of the server can be increased by controlling the arrival rate of the new incoming calls. Moreover, the system can be made more efficient by enhancing the service rate, which in turn reduces the failure frequency of the server. 8.4. Server state probabilities The efficiency of the system is also affected by the state of the sever at any instant t. Tables 3 and 4, display the variation in the server state probabilities with varying values of service rates (l1, l2) and retrial rates (h1, h2) respectively. The long run probabilities PB(t), P f ðtÞ, PR(t) and PRE2(t) increases with the growth of time t whereas PR3(t) and PRE1(t) exhibit decreasing pattern. The long run probability of the server being in busy state i.e. PB(t) increases with the rise in service rate and retrial rates. This is so because an increase in service rates makes the server more available to serve and increases its busy behavior. An increase in retrial rates for both types of calls also increases PB(t). The carried load CL(t) of the server increases with time t but decreases (increases) with the growth of service rate (retrial rate) which is very obvious. An increase in the service rate fastens the speed of the channel to handle the calls and thus reduces the load. On the contrary, carried load increases when calls kept in buffer for retry for the service in addition to new incoming calls (see Table 5). 9. Cost optimization The present section deals with the optimization of total cost incurred on the system while serving the customers. For this purpose, the cost function is framed in the following manner as:

TCðL; K; q1 ; q2 Þ ¼ CB P B ðtÞ þ Ch E½NðtÞ þ Ca PR ðtÞ þ Cb PR3 ðtÞ þ h1 C c þ h2 C d

ð66Þ

Table 4 Effect of (l1, l2) on server state probabilities. (l1, l2)

t

PB(t)

P f ðtÞ

PR(t)

PR3(t)

PRE1(t)

PRE2(t)

CL(t)

(6, 7)

1 3 5

0.5182 0.5844 0.5906

0.1872 0.3095 0.3171

0.1243 0.2768 0.2882

0.0721 0.0398 0.0351

0.0182 0.0003 0.0000

0.0541 0.0986 0.0861

9.8663 9.9399 10.0728

(6.5, 7)

1 3 5

0.5323 0.5723 0.5778

0.2135 0.3049 0.3124

0.1469 0.2698 0.2806

0.0863 0.0433 0.0391

0.0209 0.0004 0.0000

0.2003 0.1140 0.1025

8.5271 9.7725 9.8958

(8, 7)

1 3 5

0.5021 0.5408 0.5436

0.2000 0.2915 0.2982

0.1392 0.2516 0.2605

0.0767 0.0512 0.0483

0.0253 0.0008 0.0000

0.1397 0.1555 0.1476

8.7624 9.3175 9.4084

(6, 8)

1 3 5

0.4854 0.5761 0.5834

0.2266 0.3047 0.3130

0.1665 0.2725 0.2844

0.0800 0.0393 0.0347

0.0203 0.0004 0.0000

0.2339 0.1116 0.0975

8.3443 9.8798 10.0209

(6, 9)

1 3 5

0.5522 0.5686 0.5763

0.2086 0.3002 0.3091

0.1469 0.2684 0.2808

0.0775 0.0388 0.0343

0.0214 0.0004 0.0000

0.1857 0.1237 0.1085

8.7481 9.8252 9.9703

340

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

Table 5 Effect of (h1, h2) on server state probabilities. (h1, h2)

t

PB(t)

P f ðtÞ

PR(t)

PR3(t)

PRE1(t)

PRE2(t)

CL(t)

(1, 0.5)

1 3 5

0.5162 0.5341 0.5375

0.2016 0.2878 0.2947

0.1365 0.2484 0.2574

0.0851 0.0506 0.0478

0.0216 0.0005 0.0000

0.2250 0.1664 0.1574

8.2811 9.2746 9.3687

(1.5, 0.5)

1 3 5

0.5188 0.5339 0.5374

0.2021 0.2880 0.2947

0.1369 0.2484 0.2574

0.0852 0.0507 0.0478

0.0188 0.0003 0.0000

0.2243 0.1666 0.1574

8.2951 9.2739 9.3686

(2, 0.5)

1 3 5

0.5203 0.5338 0.5374

0.2025 0.2880 0.2947

0.1371 0.2485 0.2574

0.0854 0.0507 0.0478

0.0164 0.0002 0.0000

0.2243 0.1667 0.1574

8.3019 9.2733 9.3685

(2, 2)

1 3 5

0.5897 0.5509 0.5538

0.2066 0.2984 0.3040

0.1421 0.2574 0.2656

0.0827 0.0525 0.0492

0.0170 0.0004 0.0000

0.1473 0.1385 0.1314

8.6347 9.3855 9.4742

(2, 2.5)

1 3 5

0.5929 0.5553 0.5581

0.2082 0.3011 0.3065

0.1429 0.2598 0.2678

0.0839 0.0529 0.0495

0.0174 0.0004 0.0000

0.1411 0.1311 0.1246

8.6399 9.4148 9.5020

where, CB: Cost per unit time when the channel is busy; Ch: Holding cost per unit time of each call present in the system; Ca: Repair cost incurred per unit time for a broken down channel; Cb: Cost incurred per unit time for the channel being in broken down state but the repair is not yet started; Cc: Fixed cost incurred when handoff calls (priority customer) retry for the service each time; Cd: Fixed cost incurred when new calls (non-priority customer) retry for the service each time; In order to obtain the optimal values of capacity and threshold recovery parameters (L, K, q1, q2), we minimize the total cost of the system. The non-linear optimization problem is solved by using direct search approach based on discrete allocation. The optimization problem is framed mathematically as-

  ðOPÞ : TC L ; q1 ; K  ; q2 ¼ Minimize TCðL; q1 ; K; q2 Þ subject to : 1  q1  L  1 and; 1  q2  K  1 For numerical computations the default values of different cost parameters are taken as CB = 20, Ch = 15, Ca = 25, Cb = 10, Cc = 40, Cd = 40. The values for other parameters are considered as: k1 ¼ 5; k2 ¼ 6; h1 ¼ 0:5; h2 ¼ 1; l1 ¼ 8; l2 ¼ 7; b1 ¼ 1; b2 ¼ 2; a1 ¼ 0:5; a2 ¼ 0:5. Also, we fix t as 2 units and determine the optimal values of (L, K, q1, q2) using direct search approach. 9.1. Determination of optimal set (q1 ; q2 ) Now, we proceed to find out optimal set of threshold parameters (q1 ; q2 ) that produces minimum cost for the concerned queueing system. Illustration 2: Let us consider the double orbit retrial queueing model with following features: (a) The system can accommodate a maximum of 6 (i.e. L = 6) handoff calls and 8 (i.e. K = 8) new calls. (b) The threshold parameters are q1 and q2 such that (1 6 q1 6 L  1); (1 6 q2 6 K  1). (c) The values of other default parameters are fixed as: k1 ¼ 5; k2 ¼ 6; h1 ¼ 0:5; h2 ¼ 1; l1 ¼ 6; 1; b2 ¼ 2; a1 ¼ 0:5; a2 ¼ 0:5:

l2 ¼ 7; b1 ¼

In order to find the optimal values; we use ‘‘Direct Search Approach’’. In Table 6a, we vary q1 and q2 within permissible limits and search for that optimal set of threshold parameters (q1 ; q2 ) that produces minimum cost for the queueing system. Table 6a displays the result for this illustration. The optimal values for (q1 ; q2 ) are attained at as (2, 2) as shown by squared brackets with their corresponding optimal cost, average queue length and throughput. In the similar manner, different sets of optimal parameters (q1 ; q2 ) are obtained for various sets of parameters. Tables 6b– 6d display the optimal parameters obtained for different values of arrival rates (k1, k2), retrial rates (h1, h2) and repair rates (b1, b2), respectively. Conclusion: The determination of optimal threshold recovery parameters (q1 ; q2 ) can be of interest in order to initiate the start of the repair of broken down server. On this basis, cellular mobile network can handle the breakdown of channel and its repair by choosing optimal threshold parameters which may prove economical in terms of both time and money.

341

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344 Table 6a Effect of arrival rates (k1, k2) on optimal parameters (q1 ; q2 ). q1

q2

TC(t)

E[N(t)]

TP(t)

q1

q2

TC(t)

E[N(t)]

TP(t)

1 1 1 1 1 1 1

1 2 3 4 5 6 7

269.78 253.55 253.95 254.28 254.48 254.56 254.85

11.935 11.121 11.141 11.158 11.167 11.17 11.185

5.2523 5.0683 5.0843 5.0996 5.1081 5.1078 5.1123

3 3 3 3 4 4 4

4 5 6 7 1 2 3

226.88 227.36 227.33 227.94 242.07 227.7 227.62

9.8972 9.9351 9.9319 9.969 10.735 9.9826 9.9765

3.6198 3.572 3.5622 3.573 3.7173 3.4739 3.4661

2 2 2 2 2 2 2

1 2 3 4 5 6 7

239.71 224.94 225.22 225.52 225.78 225.69 226.32

10.503 9.7637 9.7811 9.7996 9.8159 9.8087 9.8476

3.9798 3.7011 3.6982 3.6912 3.6814 3.6725 3.6756

4 4 4 4 5 5 5

4 5 6 7 1 2 3

228.31 228.68 228.66 229.46 242.6 228.62 225.54

10.02 10.045 10.044 10.087 10.812 10.076 9.8734

3.4594 3.4382 3.4248 3.4775 3.5632 3.3167 3.4303

3 3 3

1 2 3

257.31 226.45 225.7

11.782 9.8784 9.823

2.6954 3.5984 3.6143

5 5 5 5

4 5 6 7

227.66 229.62 229.57 230.2

10.009 10.14 10.136 10.18

3.4261 3.2812 3.2689 3.2716

Table 6b Effect of arrival rates (k1, k2) on optimal parameters (q1 ; q2 ). (k1, k2)

(4, 5)

(4, 6)

(5, 5)

(5, 6)

(5, 7)

(6, 5)

(6, 6)

(6, 7)

q1 q2 TC(q1 ; q2 ) E[N(t)] TP(t)

4 4 184.86 7.4692 1.1814

5 2 180.02 7.231 5

4 6 153.25 5.46 0.96

2 2 224.94 9.76 3.70

5 2 220.95 9.65 2.78

2 2 231.49 10.15 3.79

3 7 171.75 6.20 3.56

3 3 222.79 9.55 3.79

Table 6c Effect of retrial rates (h1, h2) on optimal parameters (q1 ; q2 ). (h1, h2)

(0.5, 0.5)

(0.5, 1.0)

(0.5, 1.5)

(1.0, 0.5)

(1.0, 1.0)

(1.0, 1.5)

(1.5, 0.5)

(1.5, 1.0)

q1 q2

2 1 218.93 10.46 3.93

2 2 224.94 9.76 3.70

5 4 175.49 5.74 0.15

2 1 218.69 9.13 3.71

2 2 244.92 9.76 3.70

5 4 196.66 5.8192 0.17534

3 1 228.94 8.51 3.60

2 2 264.85 9.7574 3.7029

TC(q1 ; q2 ) E[N(t)] TP(t)

Table 6d Effect of repair rates (b1, b2) on optimal parameters (q1 ; q2 ). (b1, b2)

(1, 2)

(2, 2)

(3, 2)

(4, 2)

(1, 3)

(1, 4)

(2, 3)

(2, 4)

q1 q2 TC(q1 ; q2 ) E[N(t)] TP(t)

2 2 224.94 9.76 3.70

2 2 220.21 9.4885 4.1921

2 1 185.76 6.9324 5.1593

2 1 196.1 7.6528 5.1652

4 6 214.62 9.1293 3.3621

4 3 223.95 9.6608 4.5583

3 1 88.737 0.34043 6.3451

5 1 172.11 5.9801 5.4462

Illustration 3: Consider the retrial queueing system with double orbits with fixed capacity of new calls as 8 units and threshold for new calls as 4 units with the set of default parameters k1 ¼ 5; k2 ¼ 6; h1 ¼ 0:5; h2 ¼ 1; l1 ¼ 6; l2 ¼ 7; b1 ¼ 1; b2 ¼ 2; a1 ¼ 0:5; a2 ¼ 0:5: Moreover, we vary the capacity of handover calls (L) from one base station to another base station from 2 to 6 and the corresponding threshold recovery parameter (q1) from 1 to 5 (i.e. L  1). We proceed to find out the best optimal value for q1 that minimizes the cost of the system. Fig. 7 represents the variation in the optimal cost of the system with different values of ‘L’ by ‘q1’at constant t = 2 units. It is clear from the figure that a change point i.e. dip is visible corresponding to q1 = 2. However, the cost increases with increment in L.

342

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

Fig. 7. Determination of optimal parameter (q1 ; L ).

9.2. Determination of optimal parameter set (L⁄, q1 , K⁄, q2 ) and TC(L⁄, q1 , K⁄, q2 ) In the previous subsection, we have computed the optimal parameter set (q1 ; q2 ) and further corresponding optimal cost of the queueing model by keeping L and K as constant. Now, we determine optimal parameter set (L⁄, q1 , K⁄, q2 ) corresponding to the minimum cost using direct search approach by varying L, q1, K and q2 within the assumed bounds. The various bounds are as: (a) (b) (c) (d)

Capacity of the handover calls (L): 2  L  6 Capacity of the new calls (K): 2  K  8 Threshold to start the repair of channel for handover calls (q1): 1  q1  L  1 Threshold to start the repair of channel for new calls (q2): 1  q2  K  1

Different sets of optimal parameters (L⁄, q1 , K⁄, q2 ) are obtained corresponding to various sets of default parameters. We fixed l1 ¼ 6and8 for Tables 7a and 7b respectively. The set of other parameters taken are: k1 ¼ 5; k2 ¼ 6; h1 ¼ 0:5; h2 ¼ 1; l2 ¼ 7; b1 ¼ 1; b2 ¼ 2; a1 ¼ 0:5; a2 ¼ 0:5: Tables 7a and 7b summarize various optimal set for (L⁄, q1 , K⁄, q2 ) obtained for different set of (k1, k2) and (h1, h2) and default parameters. After finding the optimal threshold parameters, we find out the optimal capacities L and K corresponding to the optimal set of threshold parameters (q1 ; q2 ). Now we consider the following illustration: Illustration 4: From Table 6c, we consider the case when (h1, h2) = (1.0, 1.5) and corresponding optimal threshold parameters (q1 ; q2 ) as (5, 4). Then, we want to determine the optimal combination of total capacity of handover calls (priority customers) and new incoming calls (non-priority customers). We vary the capacity size of handover calls (L) from 6 to 12 and capacity parameter of new calls i.e. K from 5 to 8 (minimum value of L and K must be greater than their corresponding threshold repair parameters i.e. q1 = 5 and q2 = 4). Then the optimal combination of (L, K) that minimizes the cost of the system is determined as shown in (Table 8). Fig. 8 shows the variation in the total cost TC(t) of the system with varying values of L and K. It is very clear from the figure as well from the data depicted in Table 8 that the minimum cost (shown by bold letters) 92 units is obtained for (L⁄, K⁄) = (8, 7).

Table 7a Effect of (k1, k2) on optimal parameters (L⁄, q1 , K⁄, q2 ). Optimal parameters

k1 = 2, k2 = 3

k1 = 4, k2 = 6

k1 = 5, k2 = 5

k1 = 3, k2 = 4

L⁄ q1 K⁄ q2 TC(L⁄, q1 , K⁄, q2 )

5 3 6 4 21.50

6 5 3 2 135.29

4 2 8 1 81.53

6 1 8 2 48.14

Table 7b Effect of (h1, h2) on optimal parameters (L⁄, q1 , K⁄, q2 ). Optimal parameters

h1 = 1.0, h2 = 1.0

h1 = 1.5, h2 = 1.0

h1 = 2, h2 = 3

h1 = 0.5, h2 = 0.5

L⁄ q1 K⁄ q2 TC(L⁄, q1 , K⁄, q2 )

2 1 2 1 85.46

3 2 3 2 175.47

3 2 3 2 276.36

4 2 3 2 120.38

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

343

Table 8 Total cost of the system corresponding to different values of L and K. LK

5

6

7

8

6 7 8 9 10 11 12

220 219 190 231 233 235 236

233 231 166 244 245 249 247

243 242 92 255 257 258 260

253 252 262 264 266 268 275

Fig. 8. Determination of optimal parameter (L⁄, K⁄).

Fig. 9. Determination of optimal parameter (l1 ;

l1 ).

Conclusion: The optimal parameters (L⁄, q1 , K⁄, q2 ) determined can be used for the design of optimal systems by setting the buffer capacity of orbit size. A cellular mobile network with optimal capacity for both types of calls and optimal threshold parameters can be framed which may prove economical to deal with dropping of calls in real time system. 9.3. Determination of optimal service rates In Sections 8.1 and 8.2 we have computed optimal cost of the system TC (L⁄, K⁄, q1 , q2 ) by using the optimal parameters. Now, we further wish to determine optimal service rates (l1 ; l2 ) so as to know the optimal service rate at which calls must be served at the minimum cost. Illustration 5: From Table 7a, we consider set of optimal parameters (L⁄, K⁄, q1 ; q2 ) = (6, 5, 3, 2) with default parameters as k1 ¼ 4; k2 ¼ 6; h1 ¼ 0:5; h2 ¼ 1; b1 ¼ 1; b2 ¼ 2; a1 ¼ 0:5; a2 ¼ 0:5 and vary our service rates l1 from 4 to 9 units and l2 from 5 to 9 units. We proceed to find out the optimal values of both service rates which give minimum cost. In Fig. 9, we display the TC(t) corresponding to different service rates and cost values. We can observe from Table 7a that the corresponding to this particular set of optimal parameters (L⁄, K⁄,  q1 ; q2 ) = (6, 5, 3, 2) we have considered the service rates 8 and 7and the corresponding cost was TC(t) = 135.29 units. But it is remarkable to observe from the Fig. 9 that with the variation in service rates minimum cost is obtained as TC(t) = 110.52 units at l1 = l2 = 7.5 units. By taking illustration, we have demonstrated that the system can be made more economical by serving the calls or customers at optimal service rates. This gives optimal decision parameters as L⁄ = 6, K⁄ = 5, q1 = 3, q2 = 2, l1 ¼ l2 ¼ 7:5 (see Table 8). 10. Conclusions The double orbit finite retrial queue with two types of customers has been investigated. The cost optimization of the system is proposed to determine the optimal parameters. The sensitivity analysis of the model has done in context of its

344

M. Jain et al. / Applied Mathematics and Computation 253 (2015) 324–344

applications to the working of cellular mobile network. The work presented in this article is novel and seems to be very useful for the construction of optimal system designs wherein priority to one kind of traffic is given in comparison to other type of traffic. The model is realistic as it can be applicable to various congestion situations encountered in telecommunication systems, hospitals, banks, manufacturing systems which involve servicing of two types of customers under certain priority rule. The concept of threshold recovery and double orbits incorporated are not artificial as in many routine congestion situations wherein repair does not start immediately and depends upon the queue length. The realistic queueing situations corresponding to double orbits can be realize even at parties, administrative offices, and hospitals etc. where different waiting rooms or slots are available for different classes of customers. The numerical simulation facilitated give insights to the industrial engineers, system designers and decision makers to choose optimal values of various parameters under unavoidable techno-economic constraints. Acknowledgments The authors are grateful to reviewers for their worthy suggestions for the improvement of the manuscript. The second author is also thankful to MHRD for providing financial grant in the form of Senior Research Fellowship (SRF) to carry out the research work. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]

T. Yang, J.G.C. Templeton, A survey on retrial queues, Queueing Syst. 2 (3) (1987) 201–233. J.R. Artalejo, A classified bibliography of research on retrial queues: Progress in 1990–1999, Top 7 (2) (1999) 187–211. J.R. Artalejo, Accessible bibliography on retrial queues, Math. Comput. Model. 30 (3–4) (1999) 1–6. J.R. Artalejo, G.I. Falin, Standard and retrial queueing systems: a comparative analysis, Revista Matemática Complutense 15 (2002) 101–129. J.R. Artalejo, Accessible bibliography on retrial queues, Math. Comput. Model. 51 (2010) 1071–1081. G.I. Falin, J.G.C. Templeton, Retrial Queues, Chapman and Hall, London, 1997. J.R. Artalejo, A.G. Corral, Retrial Queueing Systems, Springer, Berlin, 2008. H. Li, Y.Q. Zhao, A retrial queue with a constant retrial rate, server break downs and impatient customers, Stochastic Models 21 (2005) 531–550. N.P. Sherman, J.P. Kharoufeh, An M/M/1 retrial queue with unreliable server, Oper. Res. Lett. 34 (6) (2006) 697–705. X. Yang, A. Alfa, A class of multiserver queueing system with server failures, Comput. Ind. Eng. 56 (1) (2008) 33–43. K. Wang, D. Yang, Controlling arrivals for a queueing system with an unreliable server: Newton-Quasi method, Appl. Math. Comput. 213 (1) (2009) 92– 101. M. Jain, G.C. Sharma, R. Sharma, Optimal control of (N, F) policy for an unreliable server queue with multi optional phase repair and start up, Int. J. Math. Oper. Res. 4 (2) (2012) 152–174. B.D. Choi, Y. Chang, Single server retrial queues with priority calls, Math. Comput. Model. 30 (1999) 7–32. J.R. Artalejo, A.N. Dudin, V.I. Klimenok, Stationary analysis of a retrial queue with preemptive repeated attempts, Oper. Res. Lett. 28 (2001) 173–180. I. Dimitriou, A mixed priority retrial queue with negative arrivals, unreliable server and multiple vacations, Appl. Math. Model. 37 (2013) 1295–1309. D.V. Efrosinin, O.V. Semenova, An M/M/1 system with an unreliable device and a threshold recovery policy, J. Commun. Technol. Electron. 55 (2010) 1526–1531. D. Efrosinin, A. Winkler, Queueing system with a constant retrial rate, non-reliable server and threshold-based recovery, Eur. J. Oper. Res. 210 (2011) 594–605. M. Jain, A. Bhagat, Finite population retrial queueing model with threshold recovery, geometric arrivals and impatient customers, J. Inf. Oper. Manage. 3 (2012) 162–165. G.N. Purohit, M. Jain, S. Rani, M/M/1 retrial queue with constant retrial policy, unreliable server, threshold based recovery and state dependent arrival rates, Appl. Math. Sci. 6 (2012) 1837–1846. M. Jain, A. Bhagat, Retrial queue with threshold recovery, geometric arrivals and finite capacity, Adv. Intell. Soft Comput. 131 (2012) 1039–1048. V.V. Anisimov, Averaging methods for transient regimes in overloading retrial queueing systems, Math. Comput. Model. 30 (1999) 65–78. A. Krinik, D. Marcus, R. Shiett, L.P. Chu, Transient probabilities of a single server priority queueing system, J. Statistical Plan. Inference 101 (2002) 185– 190. G.M. Leonenko, A new formula for the transient solution of the Erlang queueing model, Stat. Prob. Lett. 79 (2009) 400–406.