N type queue with Markovian flow of breakdowns

Applied Mathematics and Computation 314 (2017) 154–172

Contents lists available at ScienceDirect

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

Analysis of unreliable BMAP/PH/N type queue with Markovian flow of breakdowns Chesoong Kim a,∗, V.I. Klimenok b,c, A.N. Dudin b a

Department of Industrial Engineering, Sangji University, Wonju, Kangwon 220-702, Republic of Korea Belarusian State University, 4, Nezavisimosti Ave., Minsk 220030, Belarus c RUDN University, 6 Miklukho-Maklaya st., Moscow 117198, Russia b

a r t i c l e

i n f o

Keywords: Unreliable queueing system Batch Markov Arrival Process Breakdowns Phase-type service and repair time distribution

a b s t r a c t Unreliability of components is the inherent feature of many real world systems and its account is vital for correct prediction of performance measures of the system. Multi-server queueing model considered in this paper allows to evaluate characteristics of the systems under much more general assumptions about the probabilistic distributions describing behavior of the system than models known in existing literature. We analyse the multi-server queue with infinite buffer and the Batch Markovian Arrival Process (BMAP) of customers. The servers are identical and independent of each other. Service time of a customer has phase-type (PH) distribution. Servers are subject to breakdowns and repairs. Breakdowns occurrence moments are defined by the Markovian Arrival Process (MAP). The breakdown causes a failure of any server, which is not under repair. When a server fails the repair period starts immediately. The duration of this period has PH distribution. A customer whose service is interrupted occupies an idle server, if any, and continues his/her service. If he/she does not see an idle server, the customer goes to the buffer with some probability and permanently leaves the system with the complementary probability. We derive the constructive ergodicity condition and calculate the stationary distribution and the main performance characteristics of the system. Illustrative numerical examples are presented. © 2017 Elsevier Inc. All rights reserved.

1. Introduction Queueing. theory is well recognized mathematical tool for solving the problems of design, performance evaluation, capacity planning and optimization of various real world systems and processes in which some restricted resources are shared by numerous requests arriving at the random times. Especially wide applications this theory found in the area of electrical engineering and telecommunications. Existing literature in this theory is huge. This paper is devoted to analysis of a multiserver queueing system with unreliable servers. So, we will mention related works concerning only multi-server queues with servers breakdowns. One of the first papers where such a queue is under study is [1]. In this paper queueing model of M/M/N type is considered. The system has N identical servers and an infinite buffer. Arrival flow of customers is defined by the stationary Poisson process, service time at each server has exponential distribution. Time till breakdown of any server and its recovering time also have exponential distribution. Behavior of the system is described by two-dimensional Markov chain one component of ∗

Corresponding author at: Belarusian State University, 4, Nezavisimosti Ave., Minsk 220030, Belarus. E-mail addresses: [email protected] (C. Kim), [email protected], [email protected] (V.I. Klimenok), [email protected] (A.N. Dudin).

http://dx.doi.org/10.1016/j.amc.2017.06.035 0 096-30 03/© 2017 Elsevier Inc. All rights reserved.

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

155

which describes the number of customers in the system while the second one defines the current number of non-broken servers. Ergodicity condition for this Markov chain is derived. The problem of computation of the stationary distribution of the chain is solved using the partial generating functions and reasonings of analyticity of these functions in the unit disc of the complex plane. After that, due to the wide practical importance of analysis of multi-server queueing systems with unreliable servers, their analysis was extended in many directions. One of these directions suggests that recovering times of broken servers are not independent. There exists a finite pool of repairmen and recovering of a broken server can start only at presence of at least one repairman in this pool. If all repairmen are busy, requests for repair of the servers form a special queue. As important early work in this direction, the paper [2] should be mentioned. In the paper, again queueing model of M/M/N type is considered. All assumptions about arrivals, service, breakdowns and recovering times distribution are the same as in [1]. However, repairing of the servers is implemented by a finite number of repairmen. Behavior of the system is described by two-dimensional Markov chain with components having the same meaning as in [1]. Ergodicity condition for this Markov chain is derived. The stationary probabilities of the states of the chain are computed in the vector form. They define so called matrix-geometric distribution. This approach for analysis of two-dimensional Markov chains was further developed by authors of [2] to more general and complicated classes of Markov chains. Besides the analysis of the distribution of the Markov chain, the authors of [2] analyzed the stationary distribution of waiting time in the system. Results of numerical experiments are presented. The authors show that an approach to the analysis of such a system via computer simulation suffers in modeling the system when the number of available servers may become small even during quite short periods of time. For references to more recent papers devoted to analysis of multi-server queueing systems with unreliable servers, the papers [5–8] and [9] can be recommended. Shortcoming of models analyzed in [1] and [2], as well as in overwhelming majority of papers concerning the multiserver queueing systems with unreliable servers, from the point of view of potential applications, is imposing strict assumption about distributions of random variables defining the operation of the systems. As a rule, these distributions are assumed to be exponential. This assumption drastically simplifies mathematical analysis of the model but now rarely holds true in real world systems, especially telecommunication systems. Due to significant progress in analysis of two-dimensional Markov chains achieved first of all by Neuts, see, e.g., books [3] and [4], and his followers, the mathematical background for the analysis of queues with more complicated arrival flows and service and repair processes was created. This background allows to deal with multi-server queues with BMAP or MAP arrival process as models of arrival and breakdowns processes and the PH distribution of the service and repair times. Detailed information about a BMAP and the PH type distribution will be presented in the next section. Usage of BMAP instead of stationary Poisson process (which is a very special case of BMAP) allows to take into account not only the mean arrival rate, but also the variance of inter-arrival times and possible correlation of successive inter-arrival times. Usage of the PH distribution instead of exponential one (which is a very special case of the PH) allows to take into account not only the average service rate, but also the variance of service time. Because arrival processes in many real world systems exhibit correlation while service and repair times may by highly variable, investigation of queues of BMAP/PH/N type is important from the point view of potential applications. Likely, the first work devoted to the analysis of a reliable BMAP/M/N type queue with infinite buffer is [10]. The BMAP/PH/N/N type queue, which is essential extension of Erlang’s loss model, was analyzed in [11]. The BMAP/PH/N type queue with retrials of customers was analyzed in [12]. The literature on multi-server queueing systems with unreliable servers and arrival flow different from the stationary Poisson process and (or) service time distribution different from exponential one is not extensive. In the paper [13], published 30 years after [2], results of [2] for the M/M/N type unreliable queue are extended to the M/PH/N type queue. Such an extension was not trivial from theoretical point of view and, what is more essential, from computational point of view due to high dimension of multi-dimensional Markov chain that should be investigated to analyze the queueing model. The main advantage of our paper comparing to other literature on multi-server queueing systems with unreliable servers consists of the analysis of the system under quite general assumptions about arrival processes of customers and breakdowns and distributions of service and repair times. The system, which is close to the system under consideration, has been investigated in [14]. But the model in [14] does not have a buffer and much less effective way to construct the multidimensional Markov chain is utilized there. We study the system of BMAP/PH/N type using another way to construct the multi-dimensional Markov chain. This way provides much lower dimension of the chain what has a drastic effect on feasibility of computer implementation of the algorithm for computation of steady state probabilities and performance measures of the system. The rest of the paper is organized as follows. In Section 2, the mathematical model is described. The process of the system states as a multi-dimensional Markov chain is constructed in Section 3. A generator of this chain as a block structured matrix is presented here. An ergodicity condition for the Markov chain under study is derived in Section 4. This condition is intuitively tractable. In general case, it is necessary to solve a finite system of linear algebraic equations to check this condition. In case of exponential distribution of service and recovering times and stationary Poisson process of breakdown arrival, this condition is given in explicit form. An algorithm for the calculation of the steady state distribution of the Markov chain is described in brief in this section. Formulas for computation of performance measures of the system are given in Section 5. Section 6 contains numerical examples which give some insight into behavior of the system and justify investigation of unreliable queue with arrival flows more general than the stationary Poisson. Section 7 concludes the paper.

156

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

2. The mathematical model We consider an N-server retrial queue with Batch Markovian Arrival Process (BMAP). The BMAP is defined by its underlying process ν t , t ≥ 0, which is an irreducible continuous-time Markov chain with the finite state space {0, . . . , W },  k and the matrix generating function D(z ) = ∞ k=0 Dk z , |z| ≤ 1. The batches of customers enter the system only at the epochs of the chain ν t , t ≥ 0 transitions. The matrices Dk , k ≥ 1 (nondiagonal entries of the matrix D0 ) define the intensities of the process ν t , t ≥ 0 transitions which are accompanied by generating the k-size batch of customers. The matrix D(1) is an infinitesimal generator of the process ν t , t ≥ 0. The arrival rate (fundamental rate of the BMAP) is defined as λ = θ D (1 )e where θ is the unique solution of the system θ D(1 ) = 0, θ e = 1, and the rate of batch arrivals is defined as λb = θ (−D0 )e. Here and in the sequel e is a column vector of appropriate size consisting of 1’s, 0 is a row vector of appropriate size consisting of 0’s. The coefficient of variation cvar of intervals between batch arrivals is given by cv2ar = 2λb θ (−D0 )−1 e − 1 while the coefficient of correlation ccorr of intervals between successive batch arrivals is calculated as ccorr = (λb θ (−D0 )−1 (D(1 ) − D0 )(−D0 )−1 e − 1 )/cv2ar . For more information about the BMAP, its history, properties, special cases and related research see [15] and the survey paper by Chakravarthy [16]. It is assumed that all servers are identical and independent of each other. Service time of a customer by a server has PH type distribution with an irreducible representation (β, S). It means the following. Service time is interpreted as the time until the continuous time Markov chain mt , t ≥ 0, with the state space {1, . . . , M + 1} reaches the single absorbing state M + 1. Transitions of the chain mt , t ≥ 0, within the state space {1, . . . , M} are defined by the sub-generator S while the intensities of transitions into the absorbing state are defined by the vector S 0 = −Se. At the service beginning epoch, the state of the process mt , t ≥ 0, is chosen within the state space {1, . . . , M} according to the probabilistic row vector β. It is assumed that the matrix S + S0 β is an irreducible one. The service rate is defined as μ = −(βS−1 e )−1 . For more information about the PH type distribution, see, e.g., [3]. If an arriving batch of customers meets several servers being idle, the customers occupy the corresponding number of servers. If the number of idle fault-free servers is insufficient, the rest of the batch is placed at the end of the queue in the random order. Busy servers are subject to breakdowns and repair. The flow of breakdowns is described by the MAP. The MAP is defined by the state space {0, 1, . . . , V } of underlying process ηt , t ≥ 0, and by the matrix generating function H (z ) = H0 + H1 z, |z| ≤ 1. The rate of breakdowns is h = ϑ H1 e where ϑ is the unique solution of the system ϑ H (1 ) = 0, ϑ e = 1. An arriving breakdown from this MAP goes to any busy or idle fault-free server with equal probability and causes a failure of the corresponding server. If the breakdown meets all servers being under repair, it is ignored. When a server fails the repair period starts immediately. The duration of this period has PH type distribution with irreducible representation (γ , T) where the vector γ and the matrix T have the dimension R. Denote T 0 = −T e. The repair time of a server does not depend on the repair time of other broken-down servers and the service time of customers occupying the working servers. The repair rate is defined as τ = −(γ T −1 e )−1 . A customer whose service was interrupted occupies an idle fault-free server, if any, and continues his/her service. If he/she does not see an idle fault-free server, he/she leaves the system with some probability p and goes to the head of the queue in the buffer with complementary probability 1 − p. Our goal is to analyze the stationary behavior of the described system. To this end, we should describe the dynamics of the systems in terms of some random process. Because the considered queueing model is quite complicated, this process should be multi-dimensional and there are many alternatives for constructing such a process. In any case, the state space of this process is infinite. But the difficulty of construction and further analysis of the process drastically varies depending on the choice of the components of the process. Using the idea described in [17,18], in the next section, we construct the random process having relatively small size of the state space.

3. Process of the system states It is intuitively clear that the process of the system operation can be described in terms of multi-dimensional continuous time Markov chain. Our first goal is to properly choose the components of this chain. Because we consider the multi-server system with phase-type distribution of service and repair times, it is necessary to carefully describe the service and repair processes in the system, i.e., the processes of change in service and repair phases at the working and failed servers. Consider, for example, the service process. There are two known approaches to describe the PH service processes at N working servers. The more direct and readerfriendly approach assumes that the changes in PH service phases at N identical servers are described by the process ηt = {ηt(1) , ηt(2) , . . . , ηt(N ) }, t ≥ 0, where ηt(n) is a phase of service process at the nth busy server at time t. This approach is used in the most research devoted to investigation of queues with PH service time distribution. If the number of PH service phases is equal to M, the dimension of the state space of the process ηt is MN and becomes very large when the value of N increases. This implies significant difficulties in calculating the stationary distribution of the system states due to the extreme increase in the requirements to computer memory and computation time. To avoid these difficulties, one can use the another approach presented in [17,18]. With this approach, instead of the process ηt it is considered the process ψt = {ψt(1) , ψt(2) , . . . , ψt(M ) }, t ≥ 0, where ψt(m) is the number of servers at the phase m of service. The dimension of the state

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

157

N+M−1 )! space of the process ψ t is (N! (M−1 )! and can be much less then the dimension of the state space of the process η t . For example, if N = 25 and M = 2, the above-mentioned dimensions are equal to 225 and 26, respectively. Because of obvious advantages of the second approach, we will use this approach in the construction of the Markov chain describing the operation of the queue under consideration. Let, at time t,

• • • •

it be the number of customers in the queue, it ≥ 0; nt be the total number of busy servers and servers under repair, nt = 0, N; rt be the number of busy servers, rt = 0, nt ;  (m ) rt(m ) be the number of busy servers which are in the phase m, rt(m ) = 0, N, m = 1, M, M = rt ; m=1 rt R ( j) ( j) ( j) l = nt − rt ; • lt be the number of servers under repair which are in the phase j, lt = 0, N, j = 1, R, j=1 t

• ν t and ηt be the states of the directing process of the BMAP and the MAP correspondingly, νt = 0, W , ηt = 0, V . The process of the system states is described by a regular irreducible continuous time Markov chain with state space

 = {(i, n, r, ν, η, r (1) , . . . , r (M ) , l (1) , . . . , l (R) ), i = 0, n = 0, N, r = 0, n,

ν = 0, W , η = 0, V ,

M 

r ( j ) = r,

j=1

(1 )

{(i, n, r, ν, η, r , . . . , r ν = 0, W , η = 0, V ,

l ( j ) = n − r}



j=1

(M )

M 

R 

,l

(1 )

, . . . , l (R ) ), i > 0, n = N, r = 0, N,

r ( j ) = r,

j=1

R 

l ( j ) = N − r }.

j=1

The number of vectors in the state space with i = 0 is

K0 = (W + 1 )(V + 1 )

N  n  n=0 r=0

M−1 R−1 Cr+ M−1Cn−r+R−1 ,

and, for any fixed i > 0, the number of vectors in the state space is

K = (W + 1 )(V + 1 )

N 

M−1 R−1 Cr+ M−1CN−r+R−1 .

r=0

E.g., if the number of servers N is equal to 4 and the state spaces of the BMAP of customers, the MAP of breakdowns, the PH service time and the PH repair time consist of two elements (W = 1, V = 1, M = 2, R = 2 ), then K0 = 280, K = 140. When N = 5, the value K0 increases to 504 and the value K increases to 164. We suppose that the states of the process ξ t , t ≥ 0, are enumerated in the direct lexicographic order of the components it , nt , rt , ν t , ηt and in the reverse lexicographic order of components rt(1 ) , . . . , rt(M ) and components lt(1 ) , . . . , lt(R ) .



Let S˜ =

0 S0



O , S



T˜ =

0 T0



O . T

For use in the sequel, we introduce the matrices LN−i (N, S˜), LN−i (N, T˜ ), i = 1, N, Ai (N, S ), Ai (N, T ), i = 0, N, Pi (β ), Pi (γ ), i =

0, N − 1, which describe the transition rates or probabilities of the service process {rt(1 ) , . . . , rt(M ) } and repair process

{lt(1) , . . . , lt(R) }. Such kind of matrices were presented in the papers by V. Ramaswami and D. Lucantoni, see [17,18]. For the reader convenience, we give below a brief description of the meaning of introduced matrices. • The matrix Lk (l, S˜) contains the transition rates of the process {rt(1 ) , . . . , rt(M ) } leading to the releasing of one of the l − k busy servers (here k is the number of idle servers, l is a total number of idle and busy servers). • The matrix Lk (l, T˜ ) contains the transition rates of the process {lt(1 ) , . . . , lt(R ) } leading to the recovering (the repair completion) of one of the l − k servers under repair (here k is the number of idle servers, l is a total number of idle servers and servers under repair). • The matrix Pr (β) contains the transition probabilities of the process {rt(1 ) , . . . , rt(M ) } at the epoch of start of a new service given that there are r busy servers in the system. • The matrix Pn−r (γ ) contains the transition probabilities of the process {lt(1 ) , . . . , lt(R ) } at the epoch of start of a new repair process given that there are n − r failed servers in the system. • The matrix Ar (N − n + r, S ) contains the transition rates of the process {rt(1 ) , . . . , rt(M ) } within state space of this process without the increase or decrease of the number of busy servers (here r if the number of busy servers, N − n + r is a total number of idle and busy servers. • The matrix An−r (N − r, T ) contains the transition rates of the process {lt(1 ) , . . . , lt(R ) } within state space without the increase or decrease of the number of servers under repair (here n − r is the number of servers under repair, N − r is a total number of idle servers and servers under repair).

158

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

The algorithms for calculating the just described matrices are presented, e.g., in [20]. Now we introduce the notation that will be used in constructing the generator of the chain ξ t , t ≥ 0. Some notation will be presented with brief explanations.

1. en is a column vector of size n, consisting of 1’s. 2. I (O) is an identity (zero) matrix of appropriate dimension. When needed, we will identify the dimension of this matrix with suffix. 3. diag {al , l = 1, L } is a diagonal matrix with the diagonal entries al . 4.  and  are the symbols of the Kronecker product and sum of matrices, see, e.g., [19], and δ m,n is Kronecker’s symbol. 5. Al = A  . . .  A, l ≥ 1, A0 = 1.





l−1



l

6. Al = m=0 Inm  A  Inl−m−1 , l ≥ 1, for the matrix A having n rows. ¯ = W + 1, V¯ = V + 1, a = W ¯ V¯ . 7. W   n n ! 8. Cnm = = m!(n−m )! . m 9. B (n,r ) = Ia  Pr (β )  IC R−1

n−r+R−1

.

This

matrix

contains

the

transition

probabilities

of

the

process

χt =

{ν, η, r (1) , . . . , r (M ) , l (1) , . . . , l (R) } at the epoch of start of a new service given that there are r busy servers and

n − r failed servers in the system. 10. T0(n,r ) = Ia  IC M−1  LN−n (N − r, T˜ ). This matrix contains the transition rates of the process χ t leading to the recovering r+M−1

of one of the n − r servers under repair. 11. S0(n,r ) = Ia  LN−n (N − n + r, S˜)  IC R−1 . This matrix contains the transition rates of the process χ t leading to the releasn−r+R−1

ing of one of the r busy servers. 12. C (n,r ) = D 0  H0  Ar (N − n + r, S )  An−r (N − r, T ). D0  H (1 )  AN (N, T˜ ), if r = 0, ( N,r ) 13. C = The matrices C (n,r ) , r = 0, n, n = 1, N, contain the transition D0  H0  Ar (r, S )  AN−r (N − r, T ), if r = 1, N. rates of the process χ t within state space without an increase or decrease of the number of busy servers and servers under repair. 14. H (n,r ) = IW¯  H1  IC M−1  Pn−r (γ ). This matrix contains the transition rates of the process χ t caused by a breakdown r+M−1

arrival and an increase of the number of servers under repair from n − r to n − r + 1. Since n < N, the breakdown does not affect on the number of busy servers. 15. H (N,r ) = 1r IW¯  H1  Er,r−1  PN−r (γ ). This matrix contains the transition rates of the process χ t caused by a breakdown arrival and an increase of the number of servers under repair from N − r to N − r + 1. In this case, idle servers in the system are absent, therefore the breakdown reaches one of r busy servers with probability 1r and interrupts the service of a customer. 16. Pi, j (β ) = Pi (β )Pi+1 (β ) . . . Pj−1 (β ). This matrix contains the transition probabilities of the process χ t at the epoch of start of a new service on j − i servers given that there are i busy servers in the system. 17. Dl(n,r ) = Dl  IV¯  Pr,r+min{N−n,l } (β )  IC R−1 . This matrix contains the transition rates of the process χ t leading to arrival n−r+R−1

of l-size batch of customers and an increase of the number of busy servers from r to r + min{N − n, l }. M−1 M−1 18. Er,r−1 is a matrix of order Cr+ × Cr+ whose elements are defined as follows. The rows of this matrix are marked by M−1 M−2  the vectors ( j1 , j2 , . . . , jM ), jm = 1, r, m = 1, M, M lexicographic order, and the columns m=1 jm = r, ordered in the reverse  ˜ ˜ of this matrix are marked by the vectors ( j1 , j˜2 , . . . , j˜M ), j˜m = 1, r − 1, m = 1, M, M m=1 jm = r − 1, ordered in the reverse lexicographic order. Each vector ( j1 , j2 , . . . , jM ) is compared with each vector ( j˜1 , j˜2 , . . . , j˜M ) as follows. The vector (δ1 , δ2 , . . . , δM ) where δm = jm − j˜m , m = 1, M, is calculated. Let there exists m0 such that δm > 0 and δm = 0, m = m0 . Then the entry (Er )( j

0

) of the matrix Er is set equal to jm0 . If there exists at least one m0 such that ), ( j˜ , j˜ ,..., j˜ ) = 0. Under such a construction, the entry (Er )( j , j ,..., j ), ( j˜ , j˜ ,..., j˜ ) being

1 , j2 ,..., jM ), ( j1 , j2 ,..., jM

˜ ˜

δm0 < 0, then we put (Er )( j1 , j2 ,..., jM

˜

1

2

M

1

2

M

1

2

M

divided by r defines the probability that an arriving breakdown will be directed to one of the servers that serves a customer on phase m0 . It is clear, that the sum of entries of any row of the matrix Er being divided by r is equal to 1. Lemma 1. Infinitesimal generator Q of the Markov chain ξ t , t ≥ 0, has the following block structure

⎛

F0 ⎜Q˜−1 ⎜O Q =⎜ ⎜O ⎝ ...

F1 Q0 Q−1 O .. .

F2 Q1 Q0 Q−1 .. .

F3 Q2 Q1 Q0 .. .

··· ··· ··· ··· .. .

⎞ ⎟ ⎟ ⎟. ⎟ ⎠

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

where

⎛

C0 ⎜A1 ⎜. F0 = ⎜ ⎜.. ⎝O O

D 0,1 C1 .. . O O

⎛

D 0,2 D 1,2 .. . O O

··· ··· .. . ··· ···

⎞

⎞

D0,N−1 D1,N−1 .. . CN−1 AN

⎛

159

D0,N D1,N ⎟ ⎟ .. ⎟ + H + pH¯ + 0 , . ⎟ DN−1,N ⎠ CN

⎞

D0,N+k O ⎜D1,N+k ⎟ ⎜.. ⎟ ⎟ + δ1,k (1 − p)⎜. ⎟, k ≥ 1, Fk = ⎜. ⎝.. ⎠ ⎝O ⎠ DN,N+k HN Q˜−1 =



| AN BN−1

OK ×(K0 −K )



, Qk = DN,N+k + δ1,k (1 − p)HN , k ≥ 1.

Q0 = CN + pHN + ,

Q−1 = AN BN−1 ,

Here 0 and  are diagonal matrices which guarantee the equality Qe = 0,

⎛

T0(n,0) ⎜ S0(n,1) ⎜ O An = ⎜ ⎜ . ⎝ . . O

⎛

d (n,l ) = a

n+R−1

···

O O O .. .

O .. . O 

⎞

⎟ ⎟ ⎟, n = 1, N; ⎟ ⎠

O

···

Dl(n,1) .. . O

··· .. .

H˜ 0 O .. . O O



a

n  r=0

⎛

O ⎜ .. H¯ = ⎜ . ⎝O O

O .. . O O

⎛

OaC R−1

⎜ ⎝

Bn = ⎜

n+R−1

O .. . O

··· ··· .. . ··· ···

O H˜ 1 .. . O O

|O

Hn

··· .. . ··· ··· ×aCnR+−1R

···

⎞

O

⎟ O ⎟ ⎟, n = 0, N, l ≥ 0; .. ⎠ . (n,n ) Dl

M−1 R−1 Ct+ M−1Cmin{N,n+l }−t+R−1 , n = 0, N , l ≥ 0;

t=0

Oa ⎜O ⎜. H=⎜ ⎜ .. ⎝O O

Cn = diag{C (n,r ) , r = 0, n}, n = 0, N;

S0(n,n )

Dl(n,0)

×d (n,l )

min{N,n+l }−n−1

⎛

H˜ n =

T0(n,1) S0(n,2) .. . O

OaC R−1

⎜O ⎜ Dn,n+l = ⎜. ⎝.. O

··· ··· ··· .. .

O

O O .. . O O

O O H˜ N−1 OK

 M−1 Cr+ C R−1 ×aCnM−1 M−1 n−r+R−1 +M

⎞ ⎟ ⎟ ⎟; ⎟ ⎠

, n = 0, N − 1;

⎛

⎞

O .. . O O

Hn = diag {H (n,k ) , k = 0, n}, n = 0, N − 1;

OaC R−1

O .. ⎟ . ⎟; ⎠ O HN

N+R−1

⎜ H(N,1) ⎜ . ⎝ ..

HN = ⎜

O

O

···

O

O

O .. . O

··· .. . ···

O .. .

O .. .

H (N,N )

B (n,0)

O

···

O

⎞

O .. . O

B (n,1)

··· .. . ···

O .. .

⎟ ⎟, n = 0, N − 1. ⎠

.. . O

⎞ ⎟ ⎟ ⎟; ⎠

OaC M−1

N+M−1

B (n,n )

Proof. First, we clarify the probabilistic sense of the matrices An , Dn,n+l , Hn , H˜ , H, H¯ , Bn constituting blocks of the generator Q. The matrix An describes the transition rates of the chain ξ t from the states corresponding to the total number n of busy and failed servers to the states corresponding to n − 1 busy and failed servers.

160

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

The matrix Dn,n+l describes the transition rates of the chain ξ t from the states corresponding to the total number n of busy and failed servers to the states corresponding to n + l busy and failed servers. The matrices Hn and H˜ describe the transition rates of the chain ξ t that are accompanied by a breakdown arrival when there are n < N busy and failed servers in the system. In this case, the breakdown does not affect the busy servers because the interrupted customer continues his/her service on an idle server but the number of failed servers increases by one. The matrices Hn and H˜ constitute the matrix H which describes the transition rates of the chain ξ t that are accompanied by a breakdown arrival when there are idle servers in the system. In this case the number of failed servers increases by one and the number of busy servers does not change. The matrix H¯ describes the transition rates of the chain ξ t that are accompanied by a breakdown arrival when idle servers are absent in the system. In this case, the number of failed servers increases by one and the number of busy servers decreases by one. The matrix Bn describes the transition probabilities of the customers service process at a new service beginning epoch given that there are n busy and failed servers in the system, n = 0, N − 1. Using description of the matrices given above, it is easy to see that probabilistic sense of blocks of the generator Q are as follows. The block Fk describes the transition rates of the chain ξ t from the states corresponding to 0 customers in the queue (in this case the total number of busy and failed servers may be equal to 0, 1, . . . , N) to the states corresponding to k customers in the queue, k ≥ 0. If k = 0, the total number of busy and failed servers in the system may be equal to 0, 1, . . . , N, and, if k > N, the total number of busy and failed servers is equal to N. The block Q˜−1 describes the transition rates of the chain ξ t from the states corresponding to i = 1 customers in the queue and N busy and failed servers to the states corresponding to i = 0 customers in the queue and 0, 1, . . . , N busy and failed servers. The block Qk describes the transition rates of the chain ξ t from the states corresponding to i, i > 0, customers in the queue to the states corresponding to i + k customers, k ≥ −1. It is clear, that when the system is in these states, the total number of busy and failed servers is equal to N.  Remark 1. In the following, we will keep in mind that the matrix  is of the following form:

 = diag{Ia  (N,0) , Ia  (N,1) , . . . , Ia  (N,N ) } M−1 where the matrix (N, r) is a diagonal matrix of size Cr+ C R−1 defined by M−1 N−r+R−1

(N,r ) = −diag{(1 − δ0,r )L0 (r, S˜)Pr−1 (β )  ICN−r+ e R−1 R−1 + (1 − δN,r )[Pr (β )  L0 (N − r, T˜ )]e + [Ar (r, S )  AN−r (N − r, T )]e}, r = 0, N. Corollary 1. The Markov chain ξ t , t ≥ 0, belongs to the class of continuous time quasi-Toeplitz Markov chains, see [21], or M/G/1 type Markov chains, see [4]. Proof is trivial and follows from comparison of the generator Q with the form of a generator in the definition of continuous time quasi-Toeplitz Markov chains, or M/G/1 type Markov chains. In the further investigation of the Markov chain ξ t , t ≥ 0, we will use the results presented in [21]. 4. Ergodicity condition. Steady state distribution When the state space of the random process that describes behavior of the system is infinite, the mandatory step in analysis of such a process is validation of existence of the stationary mode of the system operation under the given set of the system parameters. In the case when such a random process is the Markov chain, the existence of the stationary mode of the system operation is equivalent to ergodicity of the Markov chain. Theorem 1. The necessary and sufficient condition for ergodicity of the Markov chain ξ t , t ≥ 0, is the fulfillment of the inequality

λ + (1 − p )

N N−1 N   (2 )  xr(1) H1 e < xr L0 (N − r, T )e + xr(3) L0 (r, S˜)e r=1

r=0

(1)

r=1

where

xr(1) = xr (IV¯  eC M−1

C R−1 r+M−2 N−r+R

(3 )

xr

), xr(2) = xr (eV¯ CM−1  ICN−r+ ), R−1 R−1 r+M

= xr (eV¯  IC M−1  eC R−1 r+M−1

N−r+R−1

),

and the vector x = (x0 , x1 , . . . , xN ) is the unique solution to the system of linear algebraic equations

xU = 0, xe = 1 where non-zero blocks of the block-three-diagonal matrix U = (Ur,r )r,r =0,N are defined by formulas

(2)

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

161

Ur,r+1 = [IV¯  Pr (β )  L0 (N − r, T˜ )], r = 0, N − 1, Ur,r = (1 − δr,0 )[IV¯  L0 (r, S˜)Pr−1 (β )  IC R−1

N−r+R−1

]

+ [δr,0 H (1 ) + (1 − δr,0 )H0 ]  [Ar (r, S˜)  AN−r (N − r, T˜ ] + IV¯  (N, r ) , r = 0, N, Ur,r−1 =

1 H1  Er,r−1  PN−r (γ ), r = 1, N. r

Proof. It follows from [21], that a necessary and sufficient condition for ergodicity of the chain ξ t , t ≥ 0, can be formulated in terms of blocks of the generator Q as follows

y

∞ 

( k + 1 )Qk e < 0

(3)

k=0

where the vector y is the unique solution to the system of linear algebraic equations ∞ 

y

Qk = 0, ye = 1.

(4)

k=−1

Inequality (3) and system (4) can be rewritten using explicit expressions for blocks Qk , k ≥ −1, given by Lemma 1. As the result, we get that the ergodicity condition is the fulfillment of the inequality



∞ 

y

kDN,N+k e + (2 − p)HN e + CN e + e +

k=1

∞ 



DN,N+k e < 0

(5)

k=1

where the vector y is the unique solution to the system of linear algebraic equations



y AN BN−1 + CN +  + HN +

∞ 



DN,N+k

= 0, ye = 1.

(6)

k=1

Because the matrix AN BN−1 + CN +  + HN +

( CN +  + HN +

∞ 

∞  k=1

DN,N+k is a generator, the following equality holds:

DN,N+k )e = −AN BN−1 e.

k=1

Using this in (5), we arrive to the following inequality:



∞ 

y



kDN,N+k e + (1 − p)HN e − AN BN−1 e < 0.

(7)

k=1

Let y be of the form

y = (θ  x0 ,

θ  x1 , . . . , θ  xN ).

M−1 where the unknown vectors xr are of size V¯ Cr+ C R−1 , r = 0, N . M−1 N−r+R−1 Substituting in (6) and (7) the expressions for the matrices DN,N+k , HN , CN , and the expression for the vector y presented above and taking into account the known relations

θ

∞ 

kDk e = λ,

k=0

∞ 

Dk e = 0, (H0 + H1 )e = 0,

k=0

we reduce inequality (7) to the following one:

λ + (1 − p )

N  r=1

−

N 

 1  H1  Er,r−1  PN−r (γ )]e − xr [IV¯  Pr (β )  L0 (N − r, T )]e r N−1

xr [

r=0

xr [IV¯  L0 (r, S˜)Pr−1 (β )  IC R−1

r=1

and system (6) to the following one:

N−r+R−1

]e < 0,

(8)

162

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

(1 − δr,0 )xr−1 [IV¯  Pr−1 (β )  L0 (N − r + 1, T˜ )] + (1 − δr,0 )xr [IV¯  L0 (r, S˜)Pr−1 (β )  IC R−1

N−r+R−1

]

+ xr {[δr,0 H (1 ) + (1 − δr,0 )H0 ]  (Ar (r, S˜)  AN−r (N − r, T˜ )) + IV¯  (N, r ) } + (1 − δr,N )xr+1

 1

r+1



H1  Er+1,r  PN−r−1 (γ ) = 0, r = 0, N.

(9)

Applying mixed product rule (see, e.g., [19]) to the left hand side of inequality (8) and using relations Pk (β )e = e, Pk (γ )e = e, we get from (8) inequality (1). Further, denoting the matrix of system (9) as U = (Ur,r )r,r =0,N , we obtain system (2). System (2) has the unique solution because the matrix of this system is an infinitesimal generator of the vector Quasi-Birth-andDeath process (see [3]) defining joint distribution of the number of busy servers, number of broken servers and the states of the processes defining arrivals of breakdowns, service of customers and repair times of broken servers in the situation when the system is overloaded.  Remark 2. Ergodicity condition (1) is quite transparent. The left hand side of this inequality is the total rate of customers arrival for the service. The first summand, λ, is the rate of customers arrival to the system while the second summand is the rate of customers returning to the buffer due to a breakdown occurrence in situation when the system is overloaded. The right hand side of inequality (1) is the total rate of customers departure from the service in situation when the system is overloaded (service rate). The first summand is the rate of customers departure from the servers due to breakdowns while the second summand is the rate of customers departure from the servers after the service completion. It is intuitively clear that the Markov chain describing queueing model under study is ergodic if and only if the total arrival rate is less than the maximum value of the total departure rate. Following traditions, we call the value

ρ = N−1 r=0



λ + (1 − p) Nr=1 xr(1) H1 e  xr(2) L0 (N − r, T )e + Nr=1 xr(3) L0 (r, S˜)e

as a load of the system. Then, ergodicity condition takes the form ρ < 1. Corollary 2. In the case of stationary Poisson flow of breakdowns and exponentially distributed service and repair times, the ergodicity condition (1) is reduced to the following inequality

λ + (1 − p)(1 − x0 )h < μ

N 

rxr + τ

r=1

N−1 

( N − r )xr

(10)

r=0

where the probabilities xr , r = 0, N, define the distribution of the number of busy servers under overload condition and are calculated as follows:

N  r!(τ h−1 )r r x r =  N N  , r = 0, N . −1 )l l=0 l l! (τ h

(11)

Proof. In the case of stationary Poisson flow of breakdowns and exponentially distributed service and repair times, we obtain that

V¯ = M = R = 1, H1 = h, Er,r−1 = 1, L0 (r, S˜) = rμ, L0 (N − r, T˜ ) = (N − r )τ . Besides, the vectors xr(n ) , n = 1, 2, 3, become scalars, xr , r = 0, N. Using this in (1) and (2), we get inequality (10) where the probabilities xr are calculated as the unique solution of the system of linear algebraic equations

−Nτ x0 + hx1 = 0,

(N − r + 1 )τ xr−1 − [h + (N − r )τ ]xr + hxr+1 = 0, r = 1, N − 1,

(12)

τ xN−1 + hxN = 0. N 

xr = 1.

(13)

r=0

It is easy to verify by direct substitution that (11) gives the unique solution of system (12) and (13).



Remark 3. Let probabilities yr , r = 0, N, define the distribution of the number of broken servers under overload condition. It is clear that yr = xN−r , r = 0, N. By introducing notation σ = hτ −1 , from (11) we derive the formula σr

yr = Nr! l , r = 0, N. σ l=0 l!

(14)

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

163

Thus, yr , r = 0, N, defines the stationary distribution of Erlang’s model M/M/N/N with arrival rate h and service rate τ . This completely agrees with intuitive reasonings that in the considered system operating under the overload condition (when servers never are idle) the distribution of the number of broken servers does not depend on parameters of the arrival and service process. It depends only on breakdowns and repair rates and is given by formula (14). In the sequel, we assume that ergodicity condition is fulfilled. Let pi be a row vector of stationary probabilities of the states of the chain ξ t , t ≥ 0, having the value i of the first component, i ≥ 0. To calculate the vectors pi , i ≥ 0, we use the numerically stable algorithm (see [21]) which has been elaborated for calculating the stationary distribution of the multidimensional continuous time quasi-Toeplitz Markov chain. The algorithm is based on censoring technique and consists of the following principal steps. Algorithm. 1. Calculate the matrix G as the minimal nonnegative solution of the matrix equation ∞ 

Qn Gn+1 = O.

n=−1

2. Calculate the matrix G0 using the equation ∞ 

Q˜−1 + (Q0 +

Qn Gn )G0 = O,

n=1

whence it follows that

G0 = − ( Q0 +

∞ 

Qn Gn )−1 Q˜−1 .

n=1

3. Calculate the matrices

Q¯ i,l =

⎧ ∞  ⎪ Fn Gn−1 Gn−2 . . . Gl , i = 0, l ≥ 0, ⎨Fl + n=l+1

∞  ⎪ Qn−i Gn−1 Gn−2 . . . Gl , i ≥ 1, l ≥ i, ⎩Ql−i + n=l+1

where Gi = G, i ≥ 1. 4. Calculate the matrices l using the recursive formulae

l = (Q¯ 0,l +

l−1 

i Q¯ i,l )(−Q¯ l,l )−1 , l ≥ 1.

i=1

5. Calculate the vector p0 as the unique solution to the system



p0 Q¯ 0,0 = 0, ∞  p0 (eK0 + l eK ) = 1. l=1

6. Calculate the vectors pl as follows: pl = p0 l , l ≥ 1. The algorithm is numerically stable because all recursions involved operate with non-negative matrices only. 5. Performance measures Having the stationary distribution pi , i ≥ 0, been calculated we can find a number of stationary performance measures of the system. Below we give some of them. Nontrivial performance measures will be presented with brief explanations.  • Mean number of customers in the queue Lqueue = ∞ i=1 i pi e. • Mean number of busy servers

Nbusy = p0 diag{Iˆn , n = 0, N}e +

∞ 

pi IˆN e,

i=1

where Iˆn = diag{rIaC M−1

C R−1 r+M−1 n−r+R−1

, r = 0, n}, n = 0, N.

• Mean number of customers in the system L = Lqueue + Nbusy . • Mean number of servers under repair

Nrepair = p0 diag{nI − Iˆn , n = 0, N}e +

∞  i=1

pi (NI − IˆN )e.

164

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

• Mean number of available servers

Nidle = p0 diag{(N − n )Ia n

M−1 R−1 r=0 Cr+M−1 Cn−r+R−1

, n = 0, N }e.

• Joint probability to see r busy servers, n − r servers under repair and i customers in the queue

p0 (n, r ) = p0 I0(n,r ) e, i = 0, r = 0, n, n = 0, N, pi (N, r ) = pi I (N,r ) e, i > 0, r = 0, N, where the matrix I0(n,r ) of size K0 and the matrix I(n,r) of size K are defined as follows

⎛

⎛

⎞

⎞ O r−1  M−1 R−1 Od (n−1) ×aC M−1 C R−1 M−1 (a Cl+ C )×aCr+ C R−1 M−1 n−r+R−1 r+M−1 n−r+R−1 M−1 n−l+R−1 ⎜ ⎟ (n,r ) l=0 ( n,r ) ⎠, I I0 = ⎝ IaC M−1 C R−1 =⎝ ⎠, IaC M−1 C R−1 r+M−1 n−r+R−1 r+M−1 n−r+R−1 O O 

d

(n−1 )

=a

n−1  l  l=0 k=0

r−1 

CkM−1 C R−1 + +M−1 n−k+R−1

k=0

CkM−1 C R−1 . +M−1 n−k+R−1

(n,r )

Note that the matrix I0 (I (n,r ) ) selects the part of the vector p0 (pi , i > 0) corresponding to r busy servers and n − r servers under repair. • Probability that an arbitrary arriving customer meets r busy servers, n − r servers under repair and i customers in the queue (n,r )

p(0a ) (n, r ) =



p0 I0

∞  k=1



kDk  IaC M−1

C R−1 r+M−1 n−r+R−1

 pi I (n,r )

p(i a ) (N − r, r ) =

e , i = 0, r = 0, n, n = 0, N,

λ ∞ 

k=1

(15)



kDk  IaC M−1

C R−1 r+M−1 n−r+R−1

e

λ

, r = 0, N , i ≥ 1.

(16)

The brief explanation of formulas (15) and (16) is as follows. The numerators are the rates of arrivals which meet r busy servers, n − r servers under repair and i customers in the buffer. While λ is the rate of all arrivals. The ratio of these rates gives the probabilities p(i a ) (n, r ), i ≥ 0. • Probability that an arbitrary arriving customer meets j available servers



N− j ∞   a) p(idle ( j ) = λ−1 p0 I0(N− j,r ) kDk  IaC M−1

C R−1 r+M−1 n− j−r+R−1

k=1

r=0

e, j = 0, N .

(17)

Formula (17) is derived by analogy with formula (15) and (16). • Probability that an arbitrary customer reaches a server immediately upon arrival

Pimm =

p0

N j=1





N− j (N− j,r )  j I (k − j )Dk  IV¯ CM−1 r=0 0 k=0

C R−1 r+M−1 N− j−r+R−1

λ

e .

(18)

The brief derivation of formula (18) can be done as follows. Under the natural assumption that position of a customer in an arriving batch is uniformly distributed, the rate of arriving customers, which have a luck to occupy a server immediately upon arrival, is calculated as N− j N   p0 I0(N− j,r )



j  k=0

j=1 r=0

∞ 

kDk + j

Dk

k= j+1

 IV¯ C M−1

C R−1 r+M−1 N− j−r+R−1

e.

Dividing this rate by the rate λ of the BMAP and taking into account the evident relation get (18). • Probability that an arbitrary customer will be lost is defined by formula

Ploss = p

p0

N

(N,r ) ˆ Hr e r=1 I0

+

∞

λ

where

Hˆ r = IW¯  H1  IC M−1

C R−1 r+M−1 N−r+R−1

.

i=1

pi

N

r=1

I (N,r ) Hˆ r e

∞

k= j+1 Dk e

=−

j k=0

Dk e, we

(19)

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

165

To derive formula (19), we used the following reasonings. The value

p0

N ∞ N    I0(N,r ) Hˆ r e + pi I (N,r ) Hˆ r e i=1

r=1

r=1

is the rate of breakdowns which arrive when there are no idle servers and there is at least one server being busy. Every such breakdown causes the loss of a customer with the probability p. So, the numerator in (19) is the rate of customers which leave the system forever and should be considered as lost ones. Then the ratio of this rate and the rate λ of input flow gives the probability Ploss . It is possible to derive another formula for Ploss :

Ploss = 1 −

p0

N

n=1

n

(n,r ) (n,r ) S0 e r=1 I0

+

∞

i=1

λ

pi

N

r=1

I (N,r ) S0(N,r ) e

.

(20)

In this formula, the numerator of the fraction is the rate of output flow of serviced customers while the denominator is the rate of input flow at the system. So, the fraction is the probability that an arbitrary customer will not be lost. The value of Ploss is calculated as a complimentary probability. 6. Numerical examples In this section we present the results of three computational experiments on calculation of the system performance measures. Experiment 1. The goal of the first experiment is to demonstrate how the input rate and correlation in BMAP impact on the system performance measures. Let us introduce three MAPs having the same fundamental rate equal to 1 but different coefficients of correlation. MAP1 has the coefficient of correlation ccorr = 0.2 and is defined by the matrices



−1.3526 0

D0 =



0 , D= −0.04391



1.3436 0.02446



0.009 , 0.01945

MAP2 has the coefficient of correlation ccorr = 0.4 and is defined by the matrices





−3.39823 0.00101

D0 =

0 , D= −0.11024



3.36283 0.01214



0.0354 , 0.09709

MAP3 is a stationary Poisson flow. It is defined by D0 = −1 and D1 = 1. This MAP has the coefficient of correlation ccorr = 0. Based on these MAPs, we construct three BMAPs which are defined by the matrices Dk , k = 0, 4. To build these matrices, we follow such a way: The matrix D0 is the same as D0 in the corresponding MAP and the rest of the matrices are calculated as Dk = Dqk−1 (1 − q )/(1 − q4 ), k = 1, 4, where q = 0.8. Then, for each BMAP, we normalize the matrices Dk , k = 0, 4, to get the BMAP with the desired fundamental rate λ. We assume that there are N = 6 identical servers in the system. The PH distribution of the service time is defined by

 β = 1,



0 ,S=



−20 0



20 , −20

i.e., this is the Erlangian distribution of order 2 with intensity μ = 10 and coefficient of variation cvar = 0.5. The flow of breakdowns is a partial case of a MAP which is named as Markov Modulated Poisson Process (MMPP). It is defined by the matrices



H0 =

−1.331867 0.021398



0.063422 , H1 = −0.128390



1.268445 0



0 . 0.106992

For this flow, h = 0.4, cvar = 2 and ccorr = 0.3. The repair time distribution is hyperexponential of order 2 and is defined by the vector γ = (0.4, 0.6 ) and sub-generator T = diag{−5.971281, −0.50718}. It has intensity τ = 0.8 and coefficient of variation 2. The probability p of leaving the system after a breakdown is assumed to be 0.4. The main performance measures as functions of the input rate λ under the different values of the coefficient of correlation ccorr are depicted in Figs. 1–6. As it is seen from Figs. 1 and 2, the value of L expectable increases with λ increasing. The rate of increasing grows with increasing λ. The value of L significantly depends on the coefficient of correlation in the input flow. Under the same value of input rate λ, the mean number L of customers in the system is greater for greater coefficient of correlation. Besides, the difference in the value of L for different coefficients of correlation increases with λ increasing. Figs. 3–5 show the behavior of the means Nbusy , Nidle and Nrepair depending on the input rate λ and coefficient of correlation.

166

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

Fig. 1. Mean number of customers in the system, L, as a function of input rate, λ, under different coefficients of correlation, ccorr = 0 and ccorr = 0.2.

Fig. 2. Mean number of customers in the system, L, as a function of input rate, λ, under different coefficients of correlation, ccorr = 0 and ccorr = 0.4.

It is clear that Nbusy increases and Nidle decreases with λ increasing. It is also clear that Nrepair does not depend on λ because breakdowns affect busy servers as well as idle servers. Not so evident and expectable are the following observations: • All values Nbusy , Nidle and Nrepair do not depend on the correlation in the arrival process. • Dependence of Nbusy and Nidle on λ is almost linear. Fig. 6 shows the dependence of probability Pimm on λ under different values of coefficient of correlation. We see that: • The probability Pimm decreases when λ increases. This observation is easily explained by the fact that the mean number of idle servers Nidle decreases with λ increasing. • The probability Pimm essentially depends on the correlation in the input flow in the range of intermediate values of λ. • Under the same value of input rate λ, the value Pimm is lower for greater coefficient of correlation. • The convexity of the curves changes with increase of ccorr . As it is seen from Fig. 7, the loss probability Ploss first significantly increases when λ grows and, then, stabilizes at some level. The stabilization occurs earlier for larger coefficients of correlation. In general, the loss probability Ploss is higher for larger coefficients of correlation. Experiment 2. Here we investigate an impact of the breakdowns rate h on the system performance measures under different values of λ. In this experiment, we assume that an input flow is the BMAP with coefficient of correlation ccorr = 0.2 defined in Experiment 1. Normalizing the matrices Dk , k = 0, 4, we get the BMAPs with different fundamental rates: λ = 25, 35, 45. The breakdowns flow is described by the MAP from Experiment 1. To vary the breakdown rate h, we normalize the matrices H0 è H1 by multiplying them by the appropriately chosen constant. The service and repair time distributions and the probability p of leaving the system are the same as in Experiment 1.

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

167

Fig. 3. Mean number of busy servers Nbusy as a function of λ under different coefficients of correlation.

Fig. 4. Mean number of idle servers Nidle as a function of λ under different coefficients of correlation. Table 1 System load ρ as a function of λ and h.

h=0 h = 1.2 h = 2.6 h=5 h = 5.5 h=6

λ = 25

λ = 35

λ = 45

0.416668 0.542889 0.681626 0.922633 0.975399 1.02899

0.583335 0.75453 0.942691 1.26954 1.3411 1.41378

0.750016 0.966171 1.20376 1.61646 1.70681 1.79857

Figs. 8–14 illustrate the dependence of the performance measures of the system on the value of h for different values of the input rate, λ = 25, 35, 45. Fig. 8 shows the behavior of the mean number L of customers in the system. It is seen from the figure that L is an increasing function of h. This function goes up faster for larger value of λ. Such a behavior of L becomes clear if we take a look at Figs. 9–11 below. We see that, when h increases, the mean number of idle servers decreases and the mean number of servers under repair increases. This implies the smaller probability for a customer to reach a server immediately after arrival, see Fig. 12, and to be served without returning to the queue. These facts lead to the increase of the queue length. We see that the curves in Fig. 8 terminate at the points with different values of the argument h. It is explained by the fact that, when h increases, the system load ρ tends to 1 faster for larger value of λ. So, the corresponding curves terminate when h becomes such as ergodicity condition ρ < 1 is violated. To make this more obvious, we give in Table 1 the values of ρ for several values of h and λ. When ρ is not less than 1, stationary probabilities of the system states do not exist and the value L is equal to infinity.

168

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

Fig. 5. Mean number of servers under repair Nrepair as a function of λ under different coefficients of correlation.

Fig. 6. Probability Pimm as a function of λ under different coefficients of correlation.

Fig. 7. Probability Ploss as a function of λ under different coefficients of correlation.

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

169

Fig. 8. L as a function of h under λ = 25; 35; 45.

Fig. 9. Nbusy as a function of h under λ = 25; 35; 45.

Fig. 10. Nidle as a function of h under λ = 25; 35; 45.

The same effect we observe looking at Figs. 9–13. Only here the terminating point for the curves with λ = 25 is not shown due to its low informativeness. In these figures we finished calculation in the point h = 3.2. As it can be seen from Fig. 12, the probability Pimm that a customer reaches a server immediately after arrival decreases when h increases. It is explained by the fact (see Fig. 11) that the mean number of fault-free idle servers, Nidle , decreases when h increases.

170

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

Fig. 11. Nrepair as a function of h under λ = 25; 35; 45.

Fig. 12. Pimm as a function of h under λ = 25; 35; 45.

Fig. 13. Ploss as a function of h under λ = 25; 35; 45.

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172

171

Fig. 14. Zoom of Fig. 13.

Fig. 15. The probability Ploss as a function of N under different values of p.

The loss probability Ploss is depicted in Fig. 13. It is expectable that this probability increases when h increases. The dependence on λ is less predictable. We see that in the range of small values of h the loss probability is larger for larger value of λ. In the range h > 0.7, we observe the opposite trend. Here the loss probability is smaller for larger value of λ. Such a dependence on λ can be explained as follows. According to formula (19), the probability Ploss is calculated as a ratio of the rate of breakdowns, which arrive when there are no idle servers and there is at least one server being busy (let us denote this rate as hloss ) and the input rate λ. When the rate of breakdowns h is small and customer arrival rate is relatively small there is a high probability that idle servers are available and a customer whose service is interrupted have a chance to receive service on one of these servers. When customer arrival rate λ increases, the number of idle servers decreases while the number of busy servers (target-servers) increases. This results the increase of the rate hloss under the same rate h. It occurs that the rate hloss may increase more quickly than λ increases. Then, according to (19), the increase of λ is accompanied by the increase of the probability Ploss . This explains the arrangement of the curves in Figs. 13 and 14 for relatively small h. However, for the relatively large h practically all servers either provide service or are under repair and the share of the servers under repair may be essentially larger than the share of the busy servers. In this case the number of target-servers becomes smaller, the rate hloss decreases and the numerator of (19) does not increase or increases more slowly than the denominator λ. This causes the increase of the probability Ploss when λ grows. This explains the change of the arrangement of the curves in Figs. 13 and 14 when h grows. More detailed arrangement of the curves in the range of small vales of λ can be seen in Fig. 14. Experiment 3. The goal of the this experiment is to investigate the dependence of the system load ρ on the number of servers N and the probability p that an interrupted customer leaves the system permanently. Also, we are interesting to understand how the values of N and p affect the loss probability Ploss . In this experiment, the input flows is the BMAP with coefficient of correlation ccorr = 0.2 defined in Experiment 1. Normalizing the matrices Dk , k = 0, 4, we get the BMAP with fundamental rate λ = 4. The values of the system load ρ for different values of N and p are given in Table 2.

172

C. Kim et al. / Applied Mathematics and Computation 314 (2017) 154–172 Table 2 System load ρ under different p and N.

p=0 p = 0.5 p=1

N=1

N=2

N=3

N=4

N=5

N=6

N=7

N=8

0.6154 0.5962 0.5769

0.2773 0.2656 0.2539

0.1726 0.1649 0.1571

0.1243 0.1186 0.1130

0.0969 0.0925 0.0881

0.0794 0.0758 0.0722

0.0673 0.0642 0.0612

0.0584 0.0557 0.0531

The probability Ploss as a function of N for values of probability p equal to 1 and 0.5 is depicted in Fig. 15. 7. Conclusion In this paper, we have analyzed an unreliable multi-server queue with quite general assumptions about the arrival processes of customers and breakdowns, service and repair time distributions what predefines an importance of the implemented analysis from the point of view of its potential applications. The stability condition for this system is presented in a simple algorithmic form. If an arrival flow of breakdowns is described by the stationary Poisson process and service and repair times have exponential distribution, stability condition is given by a simple formula. The expressions for the key performance measures of the system in terms of steady state distribution of multi-dimensional Markov chain describing dynamics of the system are derived. The behavior of these performance measures depending on the rates of arrival of customers and breakdowns as well as on the correlation in the arrival process is numerically analyzed. The profound effect of the correlation is highlighted. This explains an importance of the provided in this paper analysis of the system with quite complicated arrival process. The assumption that arrival flows are stationary Poisson ones simplifies the analysis greatly but leads to incorrect performance evaluation of any real life system modeled by the unreliable multi-server queue with correlated flows. Sometimes errors in performance evaluation are huge while correlation in arrival flows is a typical feature of flows in many telecommunication networks. Acknowledgments This publication was financially supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (Grant No. 2014K2A1B8048465), by Belarusian Republican Foundation for Fundamental Research (Grant No. F15KOR-001) and by the Ministry of Education and Science of the Russian Federation (the Agreement number 02.a03.21.0 0 08). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]

I.L. Mitrani, B. Avi-Itzhak, A many-server queue with server interruptions, Oper. Res. 163 (1968) 28–638. M.F. Neuts, D. Lucantoni, A Markovian queue with n servers subject to breakdowns and repair, Manag. Sci. 25 (1979) 49–861. M.F. Neuts, Matrix-geometric Solutions in Stochastic Models, The Johns Hopkins University Press, Baltimore, 1981. M.F. Neuts, Structured Stochastic Matrices of M/G/1 Type and Their Applications, Marcel Dekker, New York, 1989. C.C. Kuoa, S.H. Sheub, J.C. Ke, Z.G. Zhang, Reliability-based measures for a retrial system with mixed standby components, Appl. Math. Model. 38 (2014) 4640–4651. Y.L. Hsu, J.C. Ke, T.H. Liu, C.H. Wu, Modeling of multi-server repair problem with switching failure and reboot delay and related profit analysis, Comput. Ind. Eng. 69 (2014) 21–28. C.H. Wu, J.C. Ke, Multi-server machine repair problems under a (v, r ) synchronous single vacation policy, Appl. Math. Model. 38 (2014) 2180–2189. V.I. Klimenok, A.N. Dudin, A BMAP/PH/n queue with negative customers and partial protection of service, Commun. Stat. – Simul. Comput. 41 (2012) 1062–1082. A. Dudin, C.S. Kim, S. Dudin, O. Dudina, Priority retrial queueing model operating in random environment with varying number and reservation of servers, Appl. Math. Comput. 269 (2015) 674–690. M. Combe, Queueing Models with Dependence Structure, CWI, Amsterdam, 1994. V.I. Klimenok, C.S. Kim, D.S. Orlovsky, A.N. Dudin, Lack of invariant property of erlang BMAP/PH/n/0 model, Queueing Syst. 49 (2005) 187–213. L. Breuer, A.N. Dudin, V.I. Klimenok, A retrial BMAP/PN/n system, Queueing Syst. 40 (2002) 433–457. X. Yang, A.F. Alfa, A class of multi-server queueing systems with server failures, Comput. Ind. Eng. 56 (2009) 33–43. V.I. Klimenok, D.S. Orlovsky, C.S. Kim, The BMAP/PH/n retrial queue with Markovian flow of breakdowns, Eur. J. Oper. Res. 189 (2008) 1057–1072. D. Lucantoni, New results on the single server queue with a batch Markovian arrival process, Commun. Stat. – Stoch Models 7 (1991) 1–46. S. Chakravarthy, The batch Markovian arrival process: a review and future work, in: A. Krishnamoorthy (Ed.), Advances in Probability Theory and Stochastic Processes, Notable Publications, New Jersey, 2001, pp. 21–49. V. Ramaswami, Independent Markov processes in parallel, Comm. Stat. – Stoch. Models 1 (1985) 419–432. V. Ramaswami, D.M. Lucantoni, Algorithms for the multi-server queue with phase-type service, Comm. Stat. – Stoch. Models 1 (1985) 393–417. A. Graham, Kronecker Products and Matrix Calculus with Applications, Ellis Horwood, Cichester, 1981. O. Dudina, C.S. Kim, S. Dudin, Retrial queueing system with Markovian arrival flow and phase type service time distribution, Comput. Ind. Eng. 66 (2013) 360–373. V.I. Klimenok, A.N. Dudin, Multi-dimensional asymptotically quasi-toeplitz Markov chains and their application in queueing theory, Queueing Syst. 54 (2006) 245–259.