N queue with batch-size-dependent service time

European Journal of Operational Research 244 (2015) 227–239

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Stochastics and Statistics

Algorithm for computing the queue length distribution at various time epochs in DMAP/G(1,a,b)/1/N queue with batch-size-dependent service time Miaomiao Yu a,b,∗, Attahiru Sule Alfa b a b

School of Science, Sichuan University of Science and Engineering, 643000 Zigong, Sichuan, China Department of Electrical & Computer Engineering, University of Manitoba, R3T 5V6 Winnipeg, Manitoba, Canada

a r t i c l e

i n f o

Article history: Received 28 July 2014 Accepted 29 January 2015 Available online 7 February 2015 Keywords: Queueing Batch-size-dependent service Markovian arrival process Phase-type distribution Finite-buffer

a b s t r a c t This paper presents a discrete-time single-server finite-buffer queue with Markovian arrival process and generally distributed batch-size-dependent service time. Given that infinite service time is not commonly encountered in practical situations, we suppose that the distribution of the service time has a finite support. Recently, a similar continuous-time system with Poisson input process was discussed by Banerjee and Gupta (2012). But unfortunately, their method is hard to apply in the analysis of discrete-time case with versatile Markovian point process due to the fact that the difference equation governing the boundary state probabilities is more complex than the continuous one. If we follow their ideas, we will eventually find that some important joint queue length distributions cannot be computed and thus some key performance measures cannot be derived. In this paper, replacing the finite support renewal distribution with an appropriate phase-type distribution, the joint state probabilities at various time epochs (arbitrary, pre-arrival and departure) have been obtained by using matrix analytic method and embedded Markov chain technique. Furthermore, UL-type RG-factorization is employed in numerical computation of block-structured Markov chains with finitely-many levels. Some numerical examples are presented to demonstrate the feasibility of the proposed algorithm for several service time distributions. Moreover, the impact of the correlation factor on loss probability and mean sojourn time is also investigated. © 2015 Elsevier B.V. All rights reserved.

1. Introduction An efficient yet simple algorithm for computing the queue length distribution at various time epochs has been considered as one of the most critical issues in the area of queueing theory. Vast majority of studies on the performance evaluation of queueing systems are based on the queue length information. As an interesting and significant variant of the classic queue, the algorithmic analysis of the queue length distribution in bulk service system has been extensively studied over the past decades due to its wide applications in many fields, such as in transportation, manufacturing, computercommunication networks and telecommunication systems. Among some recent researches on this topic are those by Chakravarthy (1992, 1993), Laxmi and Gupta (1999), Gupta and Laxmi (2001), Gupta and Goswami (2002), Chaudhry and Gupta (2003), Gupta and Sikdar (2004), Chaudhry and Chang (2004), Chang and Choi (2005),

∗ Corresponding author at: School of Science, Sichuan University of Science and Engineering, 643000 Zigong, Sichuan, China. Tel.: +86 13540783127. E-mail address: [email protected] (M. Yu).

http://dx.doi.org/10.1016/j.ejor.2015.01.056 0377-2217/© 2015 Elsevier B.V. All rights reserved.

Skidar and Gupta (2005), Samanta, Chaudhry, and Gupta (2007), Alfa and He (2008) and Banik, Chaudhry, and Gupta (2008), both of whom considered some numerically feasible algorithm for the stationary queue length distribution in bulk service model. However, it should be noted here that almost all of the above literature invariably assumed that the service time of each batch of customer follows the same probability distribution, regardless of the value of the batch size. Such assumption is clearly unreasonable in some real world applications. It is often seen that the expected processing time of the batch with a larger number of customers should be longer than the one with a smaller number of customers. Therefore, assuming the service time of the batch to be dependent on the batch size will be more plausible in describing the behavior of the bulk service queue. Actually, the batch-size-dependent service queue has been studied by Neuts in the 1960s of the last century (see Neuts, 1967). But after many years following the pioneering work, this type of queueing system has been relatively neglected in literature and as far as we are aware only a handful of papers have been explicitly dealing with such service mechanism. Using matrix geometric approach, Curry and Feldman (1985) considered M/M(j L,K )/1 queue and obtained the joint distribution of the number of customers in the queue and the

228

M. Yu, A. S. Alfa / European Journal of Operational Research 244 (2015) 227–239

number of customers in a served batch. The work of Bar-Lev, Parlar, Perry, Stadje, and Schouten (2007) is an extension of the work done by Curry and Feldman, in the sense that a more versatile model is considered with service times that are generally distributed. Recently, employing two different methods, the queue-length distribution of the M/G(j a,b)/1 finite-buffer queue has been carried out by Chaudhry and Gai (2012) and Banerjee and Gupta (2012), respectively. Worth mentioning here in particular is the paper of Banerjee and Gupta. In their work, the embedded Markov chain technique was employed to compute the joint distribution of the number of customers in the queue and the size of departing batch immediately after a service completion epoch. Then, they also obtained the joint distribution of the number of customers in the queue and the number of customers in a served batch at arbitrary epoch by using the supplementary variable technique and considering the supplementary variable as remaining service time. In practice, these joint distributions have their own importance as they are needed to obtain the average number of customers in the system (not only in the queue) and give a clear indication of the efficiency of the server. However, the set of performance quantities that mentioned above were not reported in the work of Chaudhry and Gai. On the other hand, from the above literature overview, we may see that the existing research about the batch-size-dependent service queue has mainly focused on the continues-time case. Its discretetime counterpart has received less attention in the past few decades. Except a limited number of studies done by several Belgian scholars (Claeys, Steyaert, Walraevens, Laevens, & Bruneel, 2013; Claeys, Walraevens, Laevens, & Bruneel, 2010) and Indian scholars (Banerjee, Gupta, & Goswami, 2014), no work in this direction has come to our notice. Numerous theoretical and computational issues for discretetime batch-size-dependent service queues have been left unexplored. In this paper, we will consider a versatile finite-buffer discrete-time queue with service times that are generally distributed and dependent on the number of customers in a served batch. Further, we also assume that customers arrive into the system according to a discretetime Markovian arrival process (DMAP). Different from independent renewal arrival process, DMAP is a very good representation of the bursty and correlated traffic arising in telecommunication systems and is a rich class of point processes which contains many familiar arrival processes such as Bernoulli arrival process, discrete-time phasetype (PH) renewal process, and Markov modulated Bernoulli process (MMBP). Thus, from a theoretical perspective, our proposed model covers various special cases. Additionally, to avoid excessive delays due to postponing service until the service control limit is reached, a modified threshold policy is also included in our model. Specifically, after a service completion epoch, if the queue size n is less than the threshold a, then the server serves them one by one and if a ≤ n ≤ b, a group of size n enters service. The maximum capacity of the server is b. This implies that the size of the batch is not greater than b. Queueing model where service can be provided either to single customer or to a group of customers is useful in different real-world operational management situations to minimize the organization costs while keeping a high level of customer satisfaction. For the sake of notational simplicity, we denote this queue by DMAP/G(1,a,b)/1/N, where N is the capacity of the system excluding those under service. In the early stage of our study, we attempted to use the method presented in the paper of Banerjee and Gupta (2012) to analyze such discrete-time queue. But we found that some joint probability distributions of the number of customers in the queue and the number of customers in a served batch cannot be computed by this method. Because the joint probability distribution that mentioned above has a higher dimension than the stationary probability vector of the underlying Markov chain of the arrival process, it is almost impossible to find the joint distributions on some boundary states by using the normalization condition. Meanwhile, since the probability of an arrival and a departure occurring simultaneously is no longer equal to zero in discrete-time queues, the analysis and computation of the joint dis-

tribution will become further complicated under such situation. How to overcome these obstacles? The algorithmic method proposed by Alfa and Li (2001) gives us an inspiration. In the 1970s, Neuts (1975) pointed out that all discrete renewal distributions with finite support can be represented by discrete PH distributions. The above excellent property of the PH distribution was employed by Alfa and Li to study the GI/G/1 queue. Replacing the discrete renewal distribution with an appropriate PH distribution, they have shown that the discrete GI/G/1 system can be analyzed as a quasi-birth-and-death (QBD) process. It is indisputable that such method substantially reduces the difficulty of the model analysis. Even more importantly, due to the use of PH distribution in queueing analysis, most results obtained in their studies are numerically tractable, and these calculations can be relatively easily done by applying a suitable software package such as Maple, Mathematica or Matlab. Thus, in this paper, we will follow the idea of Alfa and Li (2001) to investigate the DMAP/G(1,a,b)/1/N queue. This paper is structured as follows. First, we describe notations and the preliminary description of the model in Section 2. In Section 3, we carry out a detailed analysis for determining the joint state probabilities at different time epochs. A variety of numerical examples are presented in Section 4. Finally, the paper ends with Section 5 where conclusions and future scope are discussed.

2. Model description and preliminaries The queueing model under consideration has the following features. The input process is described by a DMAP, where arrivals are governed by an L-states underlying Markov chain. Suppose that the DMAP is in state i, 1 ≤ i ≤ L, at the end of an arbitrary time slot. There is a transition from state i to state j without arrivals with probability cij , 1 ≤ i, j ≤ L, and there is a transition from state i to state j with an arrival     with probability dij , 1 ≤ i, j ≤ L. Let C = cij and D = dij be the L × L nonnegative matrices. Clearly, the L × L matrix C with elements cij governs transitions without arrivals, while the L × L matrix D with elements dij governs transitions with customer arrivals. The matrix C + D ˜ be the is the transition matrix of the underlying Markov chain. Let π   ˜, ˜ C+D =π stationary probability vector of the Markov chain, i.e., π π˜ eL = 1, where eL is a column vector of 1’s of size L. Then the total ˜ DeL . Let a and arrival rate λ∗ of the Markov chain is given by λ∗ = π b denote positive integers with a < b < N. Service is either single or in batches according to the following discipline. When the amount of available customers after a service completion is smaller than a, service is provided to a single customer. If the server finds at least ‘a’ customers waiting for service, say r, then service is provided to a group of customers of size min (r, b). The service time of a batch is arbitrarily distributed and dependent on the size of the batch. Let Sr (r = 1 or a ≤ r ≤ b) denote the service time of the batch of size r with probability mass function sr,k = Pr{Sr = k}, k = 1, 2, . . . , mr , sr,0 = 0. In real world applications, the size of a serving batch is always bounded. Thus, assuming that the service time Sr has a finite support, namely mr < ∞, is very reasonable in mathematical modeling. For further use in the sequel, we introduce the following notations. Let 0 denote a zero matrix of appropriate dimension. When needed, the dimension of zero matrix will be identified with a suffix. The notation I n stands for an identity matrix of dimension n. ⊗ and ⊕ are symbols of the Kronecker product and sum of matrices. Additionally, q the notation ep denotes the column vector of dimension p with 1 in the qth position and 0 elsewhere. Throughout this paper we promise  if i < j then ij = 0. Write sr = (sr,1 , sr,2 , . . . , sr,mr ). According to Alfa (2004), the density {sr , r = 1, a, . . . , b} can be represented as a discrete PH distribution in two ways, one as a remaining time and the other as an elapsed time. In this paper, we choose the remaining time representation. Let 0 denotes the transpose of 0. The distribution of random variable Sr may be expressed as a PH distribution with representation (βr , T r ),

M. Yu, A. S. Alfa / European Journal of Operational Research 244 (2015) 227–239

where for r = 1, a, a + 1, . . . , b,

given by



βr = sr = (sr,1 , sr,2 , . . . , sr,mr ), and T r =

0 I mr −1

0 0

 k = {(k, r, i, j)|r = 1, a, . . . , b, 1 ≤ i ≤ L, 1 ≤ j ≤ mr },

.

Also, we let T 0r = emr − T r emr . By direct calculation, we have T 0r , k = 1, 2, . . . , mr . Pr {Sr = k} = βr T k−1 r In discrete-time queueing system, the time axis is divided into equal intervals called slots and all queueing activities occur at the slot boundaries. Traditionally, there are two types of systems in the discrete-time case (see, e.g., Hunter, 1983; Takagi, 1993), one is the late arrival with delayed access (LAS-DA) and the other is the early arrival system (EAS). In this paper, we consider the late arrival system with delayed access i.e. service can be completed only in (t, t+ ) and arrivals can occur only in (t− , t). 3. Calculation of joint distribution 3.1. Joint distribution of the number of customers in the queue and the number of customers in a served batch at arbitrary epoch In order to study the joint distribution as well as to obtain information on some other performance measures of the queueing model,

⎛

P0, ¯ 0¯

⎜P ⎜ 0,0¯ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ P=⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝

where δ = Lm1 + L L states given by

b

r=a

¯ for k = 0, ¯ we denote the set of mr . By level 0,

¯ i)|1 ≤ i ≤ L}. 0¯ = {(0, Relating the state of the system at two consecutive epochs t− and (t + 1)− and using lexicographical sequence for the states, we find that the corresponding transition probability matrix P of the above Markov chain is of dimension L + (N + 1)δ and has a block-structured form. Since the pattern of matrix P will depend on the values of a, b and N, we cannot write it in a unified form. Here, in order to give the reader a good intuitive understanding, we only write the matrix P in the block-partitioned form with the assumption that N = 12, a = 4 and b = 6. In other cases, P can be written similarly with a simple probabilistic arguments. No matter what values of N, a and b we try, P is always a matrix with finite GI/M/1 type structure.

⎞

P0,0 ¯ P0,0

P0,1

P1,0

P1,1

P1,2

P2,1

P2,2

P2,3

P3,2

P3,3

P3,0 P4,0

P3,4 P4,4

P5,0 P6,0

229

P4,5 P5,5

P6,1 P7,1

P5,6 P6,6

P7,2 P8,2

P6,7 P7,7

P8,3 P9,3

P7,8 P8,8

P9,4 P10,4

P8,9 P9,9

P10,5 P11,5

P9,10 P10,10

P11,6

P10,11 P11,11

P12,6 we introduce the following Markov chain for the state of the system. Let Nq (t) be the number of customers waiting in the queue (excluding those in service) at time t− ; Y (t) be the number of customers being served at time t− ; Ja (t) be the phase of the arrival process at time t− ; Js (t) be the phase of the service time of the on-going service at time t− . Then {(Nq (t), Y (t), Ja (t), Js (t)) : t ≥ 0} forms a four-dimensional Markov chain with a state space

P12,12 In matrix P, the blocks for which the entries are zero are not written. For the general case, the explicit expressions for nonzero matrices are given in the following: P0, ¯ 0¯ = C,

P0,0¯

We consider Nq (t) as the level variable and others as auxiliary variables. By level k, for 0 ≤ k ≤ N, we denote the set of δ states



P0,0 = D ⊗ β1 0L×Lma ¯

⎞ C ⊗ T 01 ⎟ ⎜ ⎜ C ⊗ T0 ⎟ ⎜ a ⎟ ⎟ ⎜ ⎟ ⎜ .. ⎟, =⎜ ⎟ ⎜ . ⎟ ⎜ ⎟ ⎜ 0 ⎟ ⎜C ⊗T b−1 ⎠ ⎝ ⎛

⎛

¯ i)|1 ≤ i ≤ L}  = {(Nq (t), Ja (t)) = (0,

   ∪ Nq (t), Y (t), Ja (t), Js (t)     0 ≤ k ≤ N, r = 1, a, a + 1, . . . , b, = k, r, i, j  . 1 ≤ i ≤ L, 1 ≤ j ≤ mr

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ P11,12 ⎠

Pi,i−1

 · · · 0L×Lmb−1 0L×Lmb ,

C ⊗ T 0b C ⊗ T 01 β1

0Lm1 ×Lma · · · 0Lm1 ×Lmb−1

0Lm1 ×Lmb

C ⊗ T 0b β1

0Lmb ×Lma · · · 0Lmb ×Lmb−1

0Lmb ×Lmb

⎞

⎟ ⎜ ⎜ C ⊗ T 0β 0Lma ×Lma · · · 0Lma ×Lmb−1 0Lma ×Lmb ⎟ ⎟ ⎜ a 1 ⎟ ⎜ ⎟ ⎜ . . . . . ⎟, . . . . . =⎜ . ⎟ ⎜ . . . . ⎟ ⎜ ⎟ ⎜ 0 ⎟ ⎜ C ⊗ T β 0Lm ×Lma · · · 0Lm ×Lm 0 Lm ×Lm b−1 b−1 b−1 b−1 b ⎠ b−1 1 ⎝ i = 1, 2, . . . , a − 1,

230

M. Yu, A. S. Alfa / European Journal of Operational Research 244 (2015) 227–239

⎛ ⎜ ⎜ ⎜ ⎜ Pi,i = ⎜ ⎜ ⎜ ⎜ ⎝

C ⊗ T 1 + D ⊗ T 01 β1 0Lm1 ×Lma · · · 0Lm1 ×Lmb−1 0Lm1 ×Lmb D ⊗ T 0a β1

C ⊗ Ta

.. .

.. .

D ⊗ T 0b−1 β1

0Lmb−1 ×Lma

D⊗ ⎛

T 0b

C ⊗ T1

⎜ ⎜ 0Lma ×Lm1 ⎜ ⎜ .. Pi,i = ⎜ ⎜ . ⎜ ⎜0 ⎝ Lmb−1 ×Lm1 0Lmb ×Lm1

β1

⎞

⎛

⎟ · · · 0Lma ×Lmb−1 0Lma ×Lmb ⎟ ⎟ ⎟ .. .. .. ⎟, . ⎟ . . ⎟ ⎟ · · · C ⊗ T b−1 0Lmb−1 ×Lmb ⎠

0Lmb ×Lma · · · 0Lmb ×Lmb−1

Pi,i−b

Pi,i+1

0Lmb ×Lm1 0Lmb ×Lma ⎛ D ⊗ T 1 0Lm1 ×Lma ⎜ ⎜ 0Lma ×Lm1 D ⊗ T a ⎜ ⎜ .. .. =⎜ ⎜ . . ⎜ ⎜0 ⎝ Lmb−1 ×Lm1 0Lmb−1 ×Lma

···

0Lm1 ×Lmb−1

0Lm1 ×Lmb

···

0Lma ×Lmb−1

0Lma ×Lmb

.. .   · · · C + D ⊗ T b−1 ..

···

· · · 0Lma ×Lmb−1 ..

⎛

Pi,0

(0, 1)

⎛

(0, a)

0Lmb−1 ×Lmb  C + D ⊗ Tb ⎞

· · · D ⊗ T b−1

··· ··· ··· ..

.

···

(0, i) · · · C ⊗ T 01 βi · · · C ⊗ T 0a βi .. .

.. .

⎞ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎠

Pb,0

(0, 1) (0, a) ⎜ ⎜ 0Lm1 ×Lm1 0Lm1 ×Lma ⎜ ⎜ ⎜ 0Lma ×Lm1 0Lma ×Lma ⎜ =⎜ .. .. ⎜ . . ⎜ ⎜ ⎜ 0 ⎝ Lmb−1 ×Lm1 0Lmb−1 ×Lma 0Lmb ×Lm1

(0, b − 1)

· · · 0Lm1 ×Lmb−1 · · · 0Lma ×Lmb−1 ..

.

.. .

· · · 0Lmb−1 ×Lmb−1

0Lmb ×Lma · · · 0Lmb ×Lmb−1

· · · 0Lma ×Lmb−1 .. .

.. .

· · · 0Lmb−1 ×Lmb−1

⎞

(C + D) ⊗ T 01 βb ⎟ (C + D) ⊗ T 0a βb ⎟ ⎟

⎟ ⎟, ⎟ ⎟ 0 (C + D) ⊗ T b−1 βb ⎟ ⎠ .. .

0Lmb ×Lma · · · 0Lmb ×Lmb−1

D ⊗ T 0b βb i = b, . . . , N − 1.

Because the vector-valued Markov process {(Nq (t), Y (t), Ja (t), Js (t)) : t ≥ 0} is a finite irreducible regular Markov chain, the following limits (stationary probabilities) exist:





¯ π0,i 1 ≤ i ≤ L, ¯ = lim Pr Nq (t ) = 0, Ja (t ) = i , t→∞

πk,r,i,j = lim Pr {Nq (t) = k, Y (t) = r, Ja (t) = i, Js (t) = j} , t→∞

0 ≤ k ≤ N, r = 1, a, a + 1, . . . , b,

1 ≤ i ≤ L, 1 ≤ j ≤ mr .

For convenience, let us combine these probabilities into the rowvectors





, π0¯ = π0,1 ¯ , π0,2 ¯ , . . . , π0,L ¯   πk,r,i = πk,r,i,1 , πk,r,i,2 , . . . , πk,r,i,mr , 0 ≤ k ≤ N, r = 1, a, a + 1, . . . , b, 1 ≤ i ≤ L,

(0, i + 1) · · · (0, b − 1) (0, b) D ⊗ T 01 βi+1 · · · 0Lm1 ×Lmb−1 0Lm1 ×Lmb D ⊗ T 0a βi+1 · · · 0Lma ×Lmb−1 0Lma ×Lmb .. .

.. .

· · · C ⊗ T 0b−1 βi D ⊗ T 0b−1 βi+1 · · · 0Lmb−1 ×Lmb−1

···

0Lm1 ×Lma · · · 0Lm1 ×Lmb−1

0Lmb ×Lm1

D ⊗ T 0b βi+1 · · · 0Lmb ×Lmb−1 ⎞

(0, b) ⎟ C ⊗ T 01 βb ⎟ ⎟ ⎟ C ⊗ T 0a βb ⎟ ⎟ ⎟, .. ⎟ . ⎟ ⎟ 0 C ⊗ T b−1 βb ⎟ ⎠ C ⊗ T 0b βb

⎞

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟, .. ⎟ . ⎟ ⎟ 0Lmb−1 ×Lmb ⎟ ⎠

i = a, . . . , b − 1,

0Lmb ×Lmb 

⎛

C ⊗ T 0b βb

0Lmb ×Lmb

.. .

0Lmb ×Lma · · · C ⊗ T 0b βi

· · · 0Lmb−1 ×Lmb−1

0Lmb ×Lm1 0Lmb ×Lma · · · 0Lmb ×Lmb−1 (C + D) ⊗ T 0b βb ⎛ ⎞ 0Lm1 ×Lm1 0Lm1 ×Lma · · · 0Lm1 ×Lmb−1 D ⊗ T 01 βb ⎜ ⎟ ⎜ 0Lma ×Lm1 0Lma ×Lma · · · 0Lma ×Lmb−1 D ⊗ T 0a βb ⎟ ⎜ ⎟ ⎜ ⎟ .. .. .. .. .. ⎟, Pi,i−b+1 = ⎜ ⎜ ⎟ . . . . . ⎜ ⎟ ⎜ ⎟ 0 ⎝ 0Lmb−1 ×Lm1 0Lmb−1 ×Lma · · · 0Lmb−1 ×Lmb−1 D ⊗ T b−1 βb ⎠

i = 0, 1, . . . , N − 1, ⎞ (0, b − 1) (0, b) ⎟ 0Lm1 ×Lmb−1 0Lm1 ×Lmb ⎟ ⎟ ⎟ 0Lma ×Lmb−1 0Lma ×Lmb ⎟ ⎟ ⎟, .. .. ⎟ . . ⎟ ⎟ 0Lmb−1 ×Lmb−1 0Lmb−1 ×Lmb ⎟ ⎠

···

0Lm1 ×Lm1

⎜ ⎜ 0Lma ×Lm1 0Lma ×Lma ⎜ ⎜ .. .. =⎜ ⎜ . . ⎜ ⎜ ⎝ 0Lmb−1 ×Lm1 0Lmb−1 ×Lma

⎟ 0Lma ×Lmb ⎟ ⎟ ⎟ .. ⎟, ⎟ . ⎟ 0Lmb−1 ×Lmb ⎟ ⎠ D ⊗ Tb

D ⊗ T 0b βa · · · 0Lmb ×Lmb−1

⎜ ⎜ 0Lm1 ×Lm1 0Lm1 ×Lma ⎜ ⎜ ⎜ 0Lma ×Lm1 0Lma ×Lma ⎜ =⎜ .. .. ⎜ . . ⎜ ⎜ ⎜ 0 ⎝ Lmb−1 ×Lm1 0Lmb−1 ×Lma 0Lmb ×Lm1

.. .

.

(0, 1) (0, a) ⎜ ⎜ 0Lm1 ×Lm1 D ⊗ T 01 βa ⎜ ⎜ 0 ⎜ 0Lma ×Lm1 D ⊗ T a βa ⎜ =⎜ .. .. ⎜ . . ⎜ ⎜ ⎜ 0 0 ⎝ Lmb−1 ×Lm1 D ⊗ T b−1 βa 0Lmb ×Lm1



.. .

⎟ C ⊗ T 0a βb ⎟ ⎟ ⎟ .. ⎟, ⎟ . ⎟ ⎟ C ⊗ T 0b−1 βb ⎠

PN,N−b

· · · 0Lm1 ×Lmb−1 0Lm1 ×Lmb

⎛

Pa−1,0

0Lmb ×Lmb−1

0Lmb ×Lma · · · 0Lmb ×Lmb−1

0Lmb ×Lm1

.. .

.

.. .

⎞

i = b + 1, . . . , N − 1,

i = a − 1, . . . , N − 1, PN,N

· · · 0Lma ×Lmb−1

0Lmb ×Lma · · · 0Lmb ×Lmb−1

0Lmb ×Lm1

i = 0, 1, . . . , a − 2, ⎞ 0Lm1 ×Lma · · · 0Lm1 ×Lmb−1 0Lm1 ×Lmb ⎟ C ⊗ T a · · · 0Lma ×Lmb−1 0Lma ×Lmb ⎟ ⎟ ⎟ .. .. .. .. ⎟, . ⎟ . . . ⎟ 0Lmb−1 ×Lma · · · C ⊗ T b−1 0Lmb−1 ×Lmb ⎟ ⎠ 0Lmb ×Lma · · · 0Lmb ×Lmb−1 C ⊗ T b

 ⎛ C + D ⊗ T 1 0Lm1 ×Lma   ⎜ ⎜ 0Lma ×Lm1 C + D ⊗ Ta ⎜ ⎜ .. .. =⎜ ⎜ . . ⎜ ⎜ 0 0Lmb−1 ×Lma ⎝ Lmb−1 ×Lm1

⎜ ⎜ 0Lma ×Lm1 0Lma ×Lma ⎜ ⎜ .. .. =⎜ ⎜ . . ⎜ ⎜ ⎝ 0Lmb−1 ×Lm1 0Lmb−1 ×Lma

C ⊗ Tb

C ⊗ T 01 βb

0Lm1 ×Lma · · · 0Lm1 ×Lmb−1

0Lm1 ×Lm1



πk,r = πk,r,1 , πk,r,2 , . . . , πk,r,L , πk π

0 ≤ k ≤ N, r = 1, a, a + 1, . . . , b,   = π k,1 , π k,a , π k,a+1 , . . . , π k,b , 0 ≤ k ≤ N,   = π 0¯ , π 0 , π 1 , . . . , π N .

  It is well known that the vector π = π 0¯ , π 0 , π 1 , . . . , π N is the unique solution to the following system of linear algebraic equations

   π 0¯ , π 0 , π 1 , . . . , π N  P − I L+(N+1)δ = 0L+(N+1)δ , π0¯ , π0 , π1 , . . . , πN eL+(N+1)δ = 1.

(1)

M. Yu, A. S. Alfa / European Journal of Operational Research 244 (2015) 227–239

  In case the vector π 0¯ , π 0 , π 1 , . . . , π N has small dimension, the finite linear system (1) with the matrix P − I L+(N+1)δ can be solved directly. However, when the dimension is large, the direct solving of the system (1) can be time and resource consuming. Taking into consideration the special structure of matrix P − I L+(N+1)δ , an effective and numerically stable algorithm for solving this system is provided in Li (2010, chap. 2). In what follows, we use the UL-type RG  factorization to compute the vector π 0¯ , π 0 , π 1 , . . . , π N . Applying RG-factorization to matrix P − I L+(N+1)δ , we have

P − I L+(N+1)δ = (I L+(N+1)δ − RU ) × diag(U 0¯ , U 0 , U 1 , . . . , U N )(I L+(N+1)δ − GL ),  I L+(N+1)δ − RU ⎛ −R0¯ IL ⎜0 Iδ ⎜ ⎜0 0 ⎜ =⎜. .. ⎜ .. . ⎜ ⎝0 0 0 0

0 −R0 Iδ .. . 0 0

⎞ 0 ⎟ 0 ⎟ ⎟ 0 ⎟ ⎟, .. ⎟ . ⎟ −RN−1 ⎠ Iδ

0 ··· 0 ··· −R1 · · · .. .. . . ··· Iδ ··· 0

⎛

0 0 0 IL ⎜ −G ¯ Iδ 0 0 0,0 ⎜ ⎜ −G ¯ −G1,0 Iδ 0 ⎜ 1,0 ⎜ −G −G −G I 2,0 2,1 ¯ δ ⎜ 2, 0 =⎜ .. .. .. .. ⎜ ⎜ . . . . ⎜ ⎝−GN−1,0¯ −GN−1,0 −GN−1,1 −GN−1,2 −GN,0¯ −GN,0 −GN,1 −GN,2

⎞ ··· 0 0 ··· 0 0⎟ ⎟ ··· 0 0⎟ ⎟ ··· 0 0⎟ ⎟. .. ⎟ .. .. ⎟ . .⎟ . ⎟ ··· Iδ 0⎠ · · · −GN,N−1 I δ

U N = PN,N − I δ ,

GN,0¯ =

j = 0, 1, 2, . . . , N − 1,

−U −1 N PN,0¯ ,

Ri = −Pi,i+1 U −1 i+1 ,

z0¯ U 0¯ = 0L , zk U k = 0δ , 0 ≤ k ≤ N.

Since U k is invertible for all 0 ≤ k ≤ N (see Li & Cao, 2004), we deduce that zk = 0δ for 0 ≤ k ≤ N. Note that the censored Markov chain U 0¯ is positive recurrent, thus we let z0¯ be the stationary probability vector of the Markov chain U 0¯ . That is to say, z0¯ U 0¯ = 0L and z0¯ eL = 1. So, we   have Z = ηz0¯ , 0δ , 0δ , . . . , 0δ and

i = 0, 1, 2, . . . , N − 1,

⎞ 0 ⎟ 0 ⎟ ⎟ 0 ⎟ ⎟ .. ⎟ . ⎟ −RN−1 ⎠ Iδ

 where the scalar η is uniquely determined by π 0¯ eL + N k=0 π k eδ = 1.   Once we get the stationary probability vector π 0¯ , π 0 , π 1 , . . . , π N , the other distributions of interest such as the joint distribution of the number of customers in the queue and the number of customers  (0 ≤ k ≤ N, r = 1, a, . . . , b), distribution of the in a served batch, πk,r number of customer in the system (including the number of customers that the server is currently processing), πkS (0 ≤ k ≤ N + b),

i = 0, 1, 2, . . . , N − 1,

U 0¯ = P0, ¯ 0¯ − I L − R0¯ U 0 G0,0¯ ,  −1  Gi,j = U i Ri U i+1 Gi+1,j − Pi,j , i = 1, 2, . . . , N − 1, j = 0, 1, . . . , i − 1,   Ri U i+1 Gi+1,0¯ − Pi,0¯ , i = 0, 1, . . . , N − 1. = U −1 i

π = S k

One may note that, according to the above equations, U i , Gi,j and Ri can be determined recursively. Employing the above UL-type   RG-factorization, the equation π P − I L+(N+1)δ = 0L+(N+1)δ can be rewritten as







π I L+(N+1)δ − RU diag U 0¯ , U 0 , U 1 , . . . , U N I L+(N+1)δ − GL



= 0L+(N+1)δ .

  Since I L+(N+1)δ − GL is invertible, this equation can be simplified further as

    π I L+(N+1)δ − RU diag U 0¯ , U 0 , U 1 , . . . , U N = 0L+(N+1)δ .

    Let Z = π I L+(N+1)δ − RU and Z = z0¯ , z0 , z1 , . . . , zN . Then the above equation can be separated into two sub-equations as



0 ··· 0 ··· −R1 · · · .. .. . . ··· Iδ ··· 0

mr L  

πk,r,i,j = πk,r eLmr , r = 1, a, . . . , b, 0 ≤ k ≤ N,

i=1 j=1

U i = Pi,i − I δ − Ri U i+1 Gi+1,i ,



0 −R0 Iδ .. . 0 0

⎧ ⎨ π 0¯ = ηz0¯ , π0 = ηz0¯ R0¯ ,  ⎩ πk = π0 ki=1 Ri−1 , k = 1, 2, . . . , N,

 πk,r =

−1 R0¯ = −P0,0 ¯ U0 ,

Gi,0¯

−R0¯ Iδ 0 .. . 0 0

distribution of the number of customers in the queue, πkQ (0 ≤ k ≤ N) and the distribution of the number of customers undergoing service with the server, πrservice (r = 0, 1, a, . . . , b) can be easily obtained and given by

Then, it follows that

GN,j = −U −1 N PN,j ,



Finally, we have

and

(I L+(N+1)δ − GL )

Clearly, this leads to

⎛ IL ⎜0 ⎜   ⎜0 π0¯ , π0 , π1 , . . . , πN ⎜ ⎜. ⎜ .. ⎜ ⎝0 0   = ηz0¯ , 0δ , 0δ , . . . , 0δ ,

where



231





π0¯ , π0 , π1 , . . . , πN I L+(N+1  )δ − RU = Z,

Zdiag U 0¯ , U 0 , U 1 , . . . , U N = 0L+(N+1)δ .

πkQ = πrservice =

⎧ π0¯ eL , ⎪ ⎪ ⎪ ⎪ ⎪ πk−1,1 eLm1 , ⎪ ⎪ ⎪ ⎪ ⎪ k ⎪  ⎪ ⎪π ⎪ eLm1 + πk−r,r eLmr , k−1,1 ⎪ ⎪ ⎨ r=a b  ⎪ ⎪ ⎪ ⎪ πk−r,r eLmr , ⎪ ⎪ ⎪ r=a ⎪ ⎪ ⎪ ⎪ b ⎪  ⎪ ⎪ ⎪ πk−r,r eLmr , ⎪ ⎩



k = 0, 1 ≤ k ≤ a − 1, a ≤ k ≤ b, b + 1 ≤ k ≤ N + a, N + a + 1 ≤ k ≤ N + b,

r=k−N

π0¯ eL + π0 eδ , k = 0, , πk eδ , 1 ≤ k ≤ N, ⎧ r = 0, ⎪ ⎪ π 0¯ eL , ⎨ ⎪ ⎪ ⎩

N 

πk,r eLmr ,

r = 1, a, . . . , b.

k=0

This completes the evaluation of joint distribution of the number of customers in the queue and the number of customers in a served batch at arbitrary epoch as well as the corresponding marginal distributions.

232

M. Yu, A. S. Alfa / European Journal of Operational Research 244 (2015) 227–239

3.2. Joint distribution of the number of customers in the queue and the number of customers in a served batch at pre-arrival epoch

Without considering the phases of arrival and service, we further define

    Let Nq− (t), Y − (t), Ja− (t), Js− (t) = k, r, i, j as seen by an arbitrary customer just before his arrival. Define

π0¯− =



t→∞





t→∞

Applying Bayes’ theorem (see Samanta, Gupta, & Sharma, 2007), we have for 1 ≤ i ≤ L, − π0,i ¯ = Pr{the queue is empty, server is idle at arbitrary epoch

and customer arrival process is in phase i | customer arrival is about to occur}

π0,l ¯ dli

l=1 ⎤ m1 mr b  N  N        q q ⎣π ¯ + πk,1 I L ⊗ em1 + πk,r I L ⊗ emr ⎦DeL 0

= ⎡

L 

= − πk,r,i,j

π0,l ¯ dli

k=0 r=a q=1

l=1

=

π0,l ¯ dli , ∗ λ

l=1

π˜ DeL = Pr{the number of customers in the queue is k, server is busy with r customers at arbitrary epoch and the respectively | customer arrival is about to occur} L 

⎤ m1 mr b  N  N        q q ⎣π ¯ + πk,1 I L ⊗ em1 + πk,r I L ⊗ emr ⎦DeL 0

=

k=0 q=1

πk,r,l,j dli = π˜ DeL

L

l=1

r = 1, a, a + 1, . . . , b, ⎛

P0,0 ⎜ P1,0 ⎜ ⎜ 0 ⎜ ⎜ 0 ⎜ ⎜ . ⎜ . ⎜ . ⎜ ⎜ 0 ⎜ ⎜ Pa,0 ⎜ ⎜Pa+1,0 ⎜ P=⎜ . ⎜ .. ⎜ ⎜P ⎜ b−1,0 ⎜ P ⎜ b,0 ⎜ ⎜ 0 ⎜ ⎜ 0 ⎜ ⎜ .. ⎜ . ⎜ ⎝ 0 0

P0,1 P1,1 P2,1 0 .. . 0 Pa,1 Pa+1,1 .. . Pb−1,1 Pb,1 Pb+1,1 0 .. . 0 0

l=1

k=0 r=a q=1

πk,r,l,j dli , 0 ≤ k ≤ N, λ∗

1 ≤ j ≤ mr .

P0,2 P1,2 P2,2 P3,2 .. . 0 Pa,2 Pa+1,2 .. . Pb−1,2 Pb,2 Pb+1,2 Pb+2,2 .. . 0 0

P0,3 P1,3 P2,3 P3,3 .. . 0 Pa,3 Pa+1,3 .. . Pb−1,3 Pb,3 Pb+1,3 Pb+2,3 .. . 0 0

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

denotes the probability that an arrival finds no customers

0,i

mr   π0¯ DeL 1  − , π = π I ⊗ eqmr DeL . k,r λ∗ λ∗ q=1 k,r L

3.3. Joint distribution of the number of customers in the queue and the size of departing batch at departure epoch In this subsection, we are going to look at the queue length only at certain selected random points. We select these points at which a bulk service is about to be completed. We will find the joint distribution of the number of customers in the queue and the size of departing batch at service completion epoch by using the method of an embedded Markov chain. We now introduce some notations as follows:

• •

Nn+ : the number of customers in the queue immediately after the nth batch departure; Yn+ : the number of customers in the nth departing batch; Jn+ : the phase of the arrival process at the nth bulk service completion epoch.

˜ = { (v1 , v2 , v3 )| 0 ≤ v1 ≤ N, v2 = 1, a, a + 1, . . . , b, 1 ≤ v3 ≤ L} . 

πk,r,l,j dli

l=1

= ⎡

i=1 j=1

   Then Nn+ , Yn+ , Jn+ , n ≥ 1 forms a three-dimensional Markov chain embedded in the queue length process, whose state space is given by

phases of arrival and service processes are i and j,

L

π0¯− =

•

L

− πk,r,i,j ,

− and πk,r,i,j , we have

0 ≤ k ≤ N, r = 1, a, a + 1, . . . , b, 1 ≤ i ≤ L, 1 ≤ j ≤ mr .

k=0 q=1

π0¯−

mr L  

− in the queue and the server is idle, πk,r represents the probability that a customer sees the server is busy with r customers and k customers are waiting in the queue. Thus, according to the expressions for π ¯−

− πk,r,i,j = lim Pr Nq− (t) = k, Y − (t) = r, Ja− (t) = i, Js− (t) = j ,

L 

− − π0,i πk,r = ¯ ,

i=1

where



− − ¯ − π0,i 1 ≤ i ≤ L, ¯ = lim Pr Nq (t ) = 0, Ja (t ) = i ,

L 

P0,a−2 P1,a−2 P2,a−2 P3,a−2 .. . Pa−1,a−2 Pa,a−2 Pa+1,a−2 .. . Pb−1,a−2 Pb,a−2 Pb+1,a−2 Pb+2,a−2 .. . 0 0

P0,a−1 P1,a−1 P2,a−1 P3,a−1 .. . Pa−1,a−1 Pa,a−1 Pa+1,a−1 .. . Pb−1,a−1 Pb,a−1 Pb+1,a−1 Pb+2,a−1 .. . 0 0

P0,a P1,a P2,a P3,a .. . Pa−1,a Pa,a Pa+1,a .. . Pb−1,a Pb,a Pb+1,a Pb+2,a .. . 0 0

  Let P = Pi,j represent the one-step transition probability matrix of the above embedded Markov chain. From a macro perspective, by partitioning the state space into levels with respect to the number of customers in the queue at an arbitrary bulk service completion epoch, we may find that the corresponding matrix P is of dimension (N + 1)(b − a + 2)L, exhibiting the following block-structured form with the assumption that N − b − 1 > a (for the case of N − b − 1 ≤ a, matrix P can also be written in a similar way):

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

P0,N−b−1 P1,N−b−1 P2,N−b−1 P3,N−b−1 .. . Pa−1,N−b−1 Pa,N−b−1 Pa+1,N−b−1 .. . Pb−1,N−b−1 Pb,N−b−1 Pb+1,N−b−1 Pb+2,N−b−1 .. . PN−1,N−b−1 0

P0,N−b P1,N−b P2,N−b P3,N−b .. . Pa−1,N−b Pa,N−b Pa+1,N−b .. . Pb−1,N−b Pb,N−b Pb+1,N−b Pb+2,N−b .. . PN−1,N−b PN,N−b

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

P0,N−1 P1,N−1 P2,N−1 P3,N−1 .. . Pa−1,N−1 Pa,N−1 Pa+1,N−1 .. . Pb−1,N−1 Pb,N−1 Pb+1,N−1 Pb+2,N−1 .. . PN−1,N−1 PN,N−1

⎞ P0,N P1,N ⎟ ⎟ P2,N ⎟ ⎟ P3,N ⎟ ⎟ .. ⎟ ⎟ . ⎟ ⎟ Pa−1,N ⎟ ⎟ Pa,N ⎟ ⎟ Pa+1,N ⎟ ⎟ .. ⎟ , . ⎟ ⎟ Pb−1,N ⎟ ⎟ Pb,N ⎟ ⎟ ⎟ Pb+1,N ⎟ ⎟ Pb+2,N ⎟ ⎟ .. ⎟ . ⎟ ⎟ PN−1,N ⎠ PN,N

M. Yu, A. S. Alfa / European Journal of Operational Research 244 (2015) 227–239

where each Pi,j is a square matrix of order (b − a + 2)L × (b − a + 2)L and 0 denotes zero matrix of size (b − a + 2)L × (b − a + 2)L. The rows of each submatrix Pi,j are indexed as (i, 1, 1), . . ., (i, 1, L), (i, a, 1), . . ., (i, a, L), (i, a + 1, 1), . . ., (i, a + 1, L), . . ., (i, b, 1), . . ., (i, b, L) and columns are indexed as (j, 1, 1), . . ., (j, 1, L), (j, a, 1), . . ., (j, a, L), (j, a + 1, 1), . . ., (j, a + 1, L), . . ., (j, b, 1), . . ., (j, b, L) with i, j = 0, 1, . . . , N. In order to obtain the explicit expression of Pi,j , we first define several auxiliary matrices. k) k)   i, j be the matrix of order L × L whose (i, j)th element a(n,r Let A(n,r is the conditional probability that given a departure, which leaves at least one customer in the queue with the input process in phase i, the next departure occurs with the arrival process in phase j and there are n(n ≥ 0) arrivals during the service of a batch of size r(r = 1 or a ≤ r ≤ b) that consists of k time slots. Conditioning on the events that may occur in the first time k) slot, we can compute matrix A(n,r recursively by using the following relationships (k)

(k−1)

(k)

(k−1)

A0,r = CA0,r , An,r = CAn,r

1 ≤ k ≤ mr , (k−1)

+ DAn−1,r ,

r = 1, a, a + 1, . . . , b,

(2)

1 ≤ n ≤ k ≤ mr ,

r = 1, a, a + 1, . . . , b,

(3)

0) k) with A(0,r = I L and A(n,r = 0L×L , 0 < mr < k < n, r = 1, a, a + 1, . . . , b. Here, the first term on the right hand side of Eq. (3) indicates that there are no arrivals in the first time slot and n customers must arrive at the system in the next k − 1 time slots. On the other hand, if one customer has arrived in the first time slot, the second term on the right hand side of Eq. (3) means that there are n − 1 arrivals in the remaining k − 1 time slots. Unconditioning on the inter-departure time, we further define mr 

An,r =

k) A(n,r sr,k ,

n ≥ 0,

r = 1, a, a + 1, . . . , b.

(4)

$ mr      k−1) C ⊗ Tr = I L ⊗ βr ⊗ T k−2 A(1,r r 

k=2

% mr       k−1) A(0,r I L ⊗ T 0r + D ⊗ I mr ⊗ T k−1 r k=1

$ r −1      m k) C ⊗ Tr ⊗ T k−1 = I L ⊗ βr A(1,r r 

k=1

% mr    k−1  I L ⊗ T 0r + D ⊗ I mr C ⊗ Tr 

k=1

  " −1    = I L ⊗ βr I Lmr − C ⊗ T r D ⊗ I mr  −1  m #   × I Lmr − C ⊗ T r I Lmr − C ⊗ T r r   × I L ⊗ T 0r , r = 1, a, a + 1, . . . , b.

A0,r =

mr  k) A(0,r ⊗ βr T k−1 T 0r r

(6)

For 1 ≤ n ≤ mr and r = 1, a, a + 1, . . . , b, we define

Q n,r =

mr mr    (k−1) k−1)  k) CAn,r + DA(n−1,r ⊗ T k−1 A(n,r ⊗ T k−1 = . r r k=n

k=n

From Eq. (6), it follows that

−1   −1    D ⊗ I mr I Lmr − C ⊗ T r Q 1,r = I Lmr − C ⊗ T r  m   × I Lmr − C ⊗ T r r , r = 1, a, a + 1, . . . , b. Using properties of Kronecker product, for 2 ≤ n ≤ mr , we may obtain the following relationship

k=max(1,n)

  Clearly, An,r is also a matrix of size L × L whose (i, j)th element an,r i, j is the conditional probability that given a departure, which leaves at least one customer in the queue with the input process in phase i, there are n(n ≥ 0) arrivals during the service of a batch of size r(r = 1 or a ≤ r ≤ b) with the input process in phase j at the next departure. Since the service time has a PH representation, using Eqs. (2) and (3), matrices An,r can be expressed as

233

Q n,r =

mr   (k−1) k−1)  CAn,r + DA(n−1,r ⊗ T k−1 r k=n



= C ⊗ Tr

mr   (k)  An,r ⊗ T k−1 r k=n

+ (D ⊗ T r )

mr   (k)  An−1,r ⊗ T k−1 r k=n−1

  = C ⊗ T r Q n,r + (D ⊗ T r ) Q n−1,r ,

r = 1, a, a + 1, . . . , b.

k=1

m  r      k) = I L ⊗ βr A(0,r ⊗ T k−1 I L ⊗ T 0r r

Clearly, with the help of Q n−1,r , Q n,r can be expressed recursively as

−1   Q n,r = I Lmr − C ⊗ T r (D ⊗ T r ) Q n−1,r ,

k=1

m  r      k k−1 = I L ⊗ βr I L ⊗ T 0r C ⊗ Tr

n ≥ 2,

r = 1, a, a + 1, . . . , b.

k=1

m  r       k−1 k−1 I L ⊗ T 0r = I L ⊗ β r C ⊗ I mr C ⊗ Tr

Thus, for n ≥ 2,

k=1

An,r =

  !−1   = I L ⊗ βr C ⊗ I mr I Lmr − C ⊗ T r  " m #   × I Lmr − C ⊗ T r r I L ⊗ T 0r , r = 1, a, a + 1, . . . , b,

k=n

r = 1, a, a + 1, . . . , b. (5)

A1,r =

mr  k) A(1,r ⊗ βr T k−1 T 0r r k=1

m  r      k) = I L ⊗ βr A(1,r ⊗ T k−1 I L ⊗ T 0r r k=1



= I L ⊗ βr

m r    k=1

(k−1)

CA1,r

(k−1)

+ DA0,r



 ⊗

T k−1 r

  I L ⊗ T 0r

mr      k) A(n,r ⊗ βr T k−1 T 0r = I L ⊗ βr Q n,r I L ⊗ T 0r , r

(7) (k) 

 Let B(nk) be the matrix of order L × L whose (i, j)th element bn i, j is the conditional probability that given a departure which leaves the system empty with the arrival process in phase i, the next departure occurs after k time slots with the arrival process in phase j and there are n arrivals during the service time of the next batch departure. According to the model assumption, when the queue length is less than a, the service is only provided for a single customer based on the order of their arrival. Thus, under such a situation the size of the next departing batch is exactly one. This is why we drop the subscript ‘r’ in the above notation.

234

M. Yu, A. S. Alfa / European Journal of Operational Research 244 (2015) 227–239 k) Employing the definition of A(n,r and after some simple probabilis-

(k)

tic arguments, Bn can be written as follows

⎧ k−1 ⎪  ⎪ (k−q) ⎪ ⎪ B(0k) = C q−1 DA0,1 , ⎪ ⎪ ⎨ q=1 k−n ⎪  ⎪ ⎪ (k−q) ⎪ ⎪ B(nk) = C q−1 DAn,1 , ⎪ ⎩



k ≥ 2, k ≥ n + 1,

n = 1, 2, . . . , m1 .

q=1

Similarly, unconditioning on the inter-departure time, we also define a new matrix on the basis of B(nk),

⎧ k−1 ∞   ⎪ (k−q) ⎪ ⎪ B = C q−1 DA0,1 s1,k−q , ⎪ 0 ⎪ ⎪ ⎪ k=2 q=max(1,k−m1 ) ⎪ ⎨ ∞ k−n   (k−q) ⎪ ⎪ B = C q−1 DAn,1 s1,k−q , n = 1, 2, . . . , m1 , n ⎪ ⎪ ⎪ ⎪ k=n+1 q=max(1,k−m1 ) ⎪ ⎪ ⎩ Bn = 0L×L , n > m1 .

Now, we are facing the same problem. It is well-known that the calculation of higher-dimensional matrices is rather complex and sometimes even harder to implement due to the computational limitations and memory requirements. In order to efficiently solve the system of Eq. (9), we still use the UL-type RG-factorization, given in Li (2010), to find the solution of Eq. (9). The major advantage of RG-factorization is that it can avoid the calculation of high dimensional matrices by decomposing them into small ones. Here, the UL-type RG-factorization of matrix P − I (N+1)(b−a+2)L is given by

(I (N+1)(b−a+2)Lδ − Rˆ U ) ⎛

ˆ 0,1 −R ˆ 0,2 −R ˆ 0,3 ··· I (b−a+2)L −R ⎜ ˆ ˆ 0 I (b−a+2)L −R1,2 −R1,3 ··· ⎜ ⎜ ˆ 2,3 0 0 I (b−a+2)L −R ··· ⎜ =⎜ .. .. .. ⎜ .. .. ⎜ . . . . . ⎜ ⎝ 0 0 0 · · · I (b−a+2)L 0 0 0 ··· 0

(8)

⎧ ! ⊗ Bj 0(b−a+2)L×(b−a+1)L , e j = 0, 1, . . . , N − 1, ⎪ ⎪ ⎨ $ b−a+2 % m1  = ⎪ e ⊗ Bk 0(b−a+2)L×(b−a+1)L , j = N. ⎪ ⎩ b−a+2

(I (N+1)(b−a+2)L − Gˆ L ) ⎛

I (b−a+2)L ⎜ −G ˆ 1,0 ⎜ ⎜ ˆ 2,0 ⎜ −G ⎜ ˆ ⎜ = ⎜ −G3,0 ⎜ .. ⎜ . ⎜ ⎜ ˆ ⎝ −GN−1,0 ˆ N,0 −G

k=N

! ⎧ eb−a+2 ⊗ Aj−i+1,1 0(b−a+2)L×(b−a+1)L , ⎪ ⎪ ⎪ ⎪ ⎨ j + 1 ≥ i, j = 0, 1, . . . , N − 1, % Pi,j = $ m1  ⎪ ⎪ ⎪ Ak,1 0(b−a+2)L×(b−a+1)L , j = N. ⎪ ⎩ eb−a+2 ⊗ k=N−i+1

Pij =

⎪ 0(b−a+2)L×(i−a+1)L ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

For i = b, b + 1, . . . , N,

⎧ 0(b−a+2)L×(b−a+1)L ⎪ ⎪ ⎪ ⎪ ⎨ Pij = $ ⎪ ⎪ ⎪ ⎪ 0(b−a+2)L×(b−a+1)L ⎩

eb−a+2 ⊗ Aj,i

0(b−a+2)L×(b−i)L ,

j = 0, 1, . . . , N − 1, eb−a+2 ⊗



Ak,i

0(b−a+2)L×(b−i)L

k=N

j = 0, 1, . . . , N − 1,

−1 ˆN , −Pi,N U

i = 0, 1, . . . , N − 1,

ˆ i = Pi,i − I (b−a+2)L − U ⎛ ˆ i,j = G

−1 ˆi U

⎝

N 

ˆ k,i , ˆ i,k U ˆ kG R

k=i+1 N 

i = 0, 1, . . . , N − 1,

⎞

ˆ k,j − Pi,j ⎠ , ˆ i,k U ˆ kG R

k=i+1

! eb−a+2 ⊗ Aj−(i−b),b , j ≥ i − b, j = 0, 1, . . . , N, % mb  eb−a+2 ⊗ Ak,b , j = N. k=N−(i−b)





π+v1 = π+v1 ,1 , π+v1 ,a , π+v1 ,a+1 , . . . , π+v1 ,b , v1 = 0, 1, . . . , N, 

π+v1 ,v2 = πv+1 ,v2 ,1 , πv+1 ,v2 ,2 , . . . , πv+1 ,v2 ,L , v1 = 0, 1, . . . , N, v2 = 1, a, a + 1, . . . , b,

0 0 0

ˆ N,j = −U ˆ N PN,j , G

,

j = N.

0 0

ˆ N = PN,N − I (b−a+2)L , U

%

mi

0

Thus, we have

ˆ i,N = R

Since we deal with a finite and irreducible Markov chain, then it must be positive recurrent. As a result, its stationary probability vector   π+ , partitioned as π+ = π+0 , π+1 , . . . , π+N , always exists, where



⎞ ··· 0 0 ⎟ I (b−a+2)L ··· 0 0 ⎟ ⎟ ˆ 2,1 I (b−a+2)L −G ··· 0 0 ⎟ ⎟ ˆ 3,1 −G ˆ 3,2 I (b−a+2)L · · · ⎟ −G 0 0 ⎟. ⎟ .. .. .. .. .. .. ⎟ . . . . . . ⎟ ⎟ ˆ ˆ ˆ ⎠ −GN−1,1 −GN−1,2 −GN−1,3 · · · I (b−a+2)L 0 ˆ N,1 −G ˆ N,2 −G ˆ N,3 · · · −G ˆ N,N−1 I (b−a+2)L −G

−1

!

⎞ ˆ 0,N −R ˆ 1,N ⎟ −R ⎟ ˆ 2,N ⎟ −R ⎟ ⎟, .. ⎟ ⎟ . ⎟ ˆ N−1,N ⎠ −R I (b−a+2)L

and

For i = 1, 2, . . . , a − 1,

For i = a, a + 1, . . . , b − 1,

(9)

where

n = 0, 1, . . . , m1 .

⎧ 0(b−a+2)L×(i−a+1)L ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨$



ˆ L ), ˆ 1, . . . , U ˆ N )(I (N+1)(b−a+2)L − G ˆ 0, U × diag(U

With the help of Eqs. (5)–(8), the submatrices Pi,j appearing in matrix P are described in detail as below: For i = 0,

P0,j



π+ P − I (N+1)(b−a+2)L = 0(N+1)(b−a+2)L , π+ e(N+1)(b−a+2)L = 1.

ˆU ) P − I (N+1)(b−a+2)L = (I (N+1)(b−a+2)L − R

Bn (n = 0, 1, 2, . . . , m1 ) is a matrix of order L × L whose (i, j)th element is the conditional probability that given a departure which leaves the system empty with the arrival process in phase i, there are n arrivals during the service time of next departure and arrival process is in  r A(k) s , changing the order of the phase j. Since An,r = m k=max(1,n) n,r r,k above double summation, we have

 −1 Bn = I − C DAn,1 ,

and πv+1 ,v2 ,v3 = limn→∞ Pr{Nn+ = v1 , Yn+ = v2 , Jn+ = v3 }. The unknown quantities {πv+1 ,v2 ,v3 } can be obtained by solving the following system of linear algebraic equations

i = 1, 2, . . . , N − 1, 0 ≤ j ≤ i − 1, ⎞ ⎛ N  −1 ˆ k,j − Pi,j ⎠ U ˆ i,k U ˆ kG ˆ i,j = ⎝ ˆj , R R k=j+1

i = 1, 2, . . . , N − 2,

i + 1 ≤ j ≤ N − 1.

ˆ i,j , i = ˆ i , i = 0, 1, . . . , N}, {G Obviously, the matrix sequences {U ˆ 1, . . . , N, j = 0, 1, . . . , N − 1} and {Ri,j , i = 0, 1 . . . , N − 1, j = 1, . . . , N} can be obtained through the backward recursions. By means of the UL-type RG-factorization, Li (2010) has provided the expression for + + the stationary probability vector (π + 0 , π 1 , . . . , π N ) as follows:



π+0 = τ x0 ,  +ˆ π+k = k−1 i=0 π i Ri,k , 1 ≤ k ≤ N

M. Yu, A. S. Alfa / European Journal of Operational Research 244 (2015) 227–239

235

Table 1 PH representation for service time of batch size r. r

βr

1

1  e 5 5

5

1  e 12 12

Tr   01×4 0 04×1 I  4  01×11 0 I 11 011×1

E [Sr ] 3

Var [Sr ] 2

6.5

11.92

where x0 is the stationary probability vector of the censored Markov ˆ 0 to level 0 and the scalar τ is uniquely determined by chain U N + π k=0 k e(b−a+2)L = 1. Having computed the stationary probability + + vector (π + 0 , π 1 , . . . , π N ), the other distributions such as the joint distribution of the number of customers in the queue and the size of de+ parting batch at departure epoch, πk,r (0 ≤ k ≤ N, r = 1, a, . . . , b) and the distribution of the number of customers in a departing batch, πrdeparture (r = 1, a, . . . , b) can be easily written as + πk,r = π+ e , 0 ≤ k ≤ N, r = 1, a, . . . , b, k,r L

πrdeparture =

N 

Lq =

k=0

kπ

Ls =

k=0

kπ

1  e 11 11

6

1  e 13 13

Tr   01×10 0 010×1 I  10  01×12 0 I 12 012×1

E [Sr ]

Var [Sr ]

6

10

7

14

loss probability and mean sojourn time, we also compare our model with an equivalent model that does not capture the correlation property. By this comparison, we can see the impact of autocorrelation on performance measures of batch-size-dependent service queue. All the calculations are performed on a PC having Intel Core i5 processor at 2.6 gigahertz with 4 gigabytes DDR3 RAM using Matlab program in Windows 7 environment. The numerical results are reported in tables up to six decimal places because of lack of space.

C=

0.3 0.1

0.1 0.1

'

&

,

D=

0.3 0.3

0.3 0.5

'

,

and this leads to λ∗ = 0.7. The distributions of Sr , r = 1, 4, 5, 6, are given by

The basic objective of modeling of queueing problem is to evaluate some queueing measures to characterize the system behavior which helps the designer of the system in taking appropriate decisions. Thus, the main goal of our investigation in this subsection is to predict various performance metrics in terms of joint probability of the number of customers in the queue and the number of customers in a served batch. Some indices to characterize the system performance can be derived from these steady state probabilities. The expected number of customers in the queue (Lq ), the average number of customers in the system (Ls ) and the average number of customers in a served batch (Lservice ) are given by N+b 

4

&

π+k,r eL , r = 1, a, . . . , b.

3.4. Performance measures

Q , k

βr

Example 1. In this example the numerical results are given for the DMAP/G(1,4,6)/1/12 queue. We choose arrival matrices as

k=0

N 

r

S k,

Lservice = π

service 1

+

b 

rπ

service . r

r=a

A particular important measure of the system performance for a finite buffer queueing system is the loss or blocking probability (Ploss ) which  − − + br=a πN,r . Based on the loss probability mentioned is equal to πN,1 above, using Little’s rule, the average waiting time of a customer in the queue (Wq ), namely the time that elapses between his arrival and the   start of his service, is given by Wq = Lq /λ , where λ = λ∗ 1 − Ploss is the effective arrival rate. Similarly, let Ws denote the mean sojourn time of a customer in the system, we also have Ws = Ls /λ . 4. Numerical examples To demonstrate the feasibility of developed algorithm obtained in the previous section, prove the importance of the analytic results and show numerically some interesting features of the system under consideration we present the results of four numerical experiments. Firstly, steady state joint probabilities at departure, arbitrary and prearrival epochs are shown in three self explanatory tables. Then, under a special discrete-Markovian arrival process, we further illustrate experimentally that the algorithm presented in this paper is valid, accurate and very reliable. Further, the effects of buffer capacity on loss probability and average waiting time in the queue are investigated under different service time distributions. On the other hand, it is generally known that the discrete-Markovian arrival process captures both burstiness and correlation of the arrival process which are highly essential to be considered to measure the performance characteristics. Thus, to show the impact of burstiness and correlation on the

1 1 , j = 1, 2, . . . , 5, Pr {S4 = j} = , 5 11 j = 1, 2, . . . , 11, 1 1 , j = 1, 2, . . . , 12, Pr {S6 = j} = , Pr {S5 = j} = 12 13 j = 1, 2, . . . , 13. Pr {S1 = j} =

Clearly, there is a stochastic order with S1 ≤st S4 ≤st S5 ≤st S6 . Such assumption is much closer to the actual service situation. The PH representations for different service times Sr are listed in Table 1. We carry out numerical work based on the procedure discussed in Section 3. In Table 2, we have computed the joint distribution of the number of customers in the queue and the size of the departing batch at departure epoch. The last column of Table 2 displays the distribution of the number of customers in the queue at departure epoch, and the last row of Table 2 presents the distribution of the number of customers in the departing batch. Then, in Tables 3 and 4, the joint distribution of the number of customers in the queue and the number of customers in a served batch at arbitrary and pre-arrival epochs are reported. Also, the marginal distributions of the number of customers in the queue and the number of customers undergoing

Table 2 Joint distribution of the number of customers in the queue and the size of departing batch at departure epoch for DMAP/G(1,4,6)/1/12 queue. k

0 1 2 3 4 5 6 7 8 9 10 11 12 Sum

r=1

r=4

r=5

r=6

Total

0.009981 0.042790 0.079495 0.101241 0.081580 0.052598 0.023825 0.005785 0.000000 0.000000 0.000000 0.000000 0.000000 0.397295

0.005769 0.018418 0.018388 0.018339 0.018139 0.017502 0.015961 0.013104 0.009075 0.004862 0.001762 0.000320 0.000000 0.141639

0.004270 0.013686 0.013669 0.013652 0.013581 0.013331 0.012652 0.011214 0.008839 0.005818 0.002942 0.001002 0.000171 0.114827

0.002978 0.011834 0.018685 0.024099 0.028338 0.031395 0.034253 0.036652 0.034746 0.031483 0.026875 0.021466 0.043435 0.346239

0.022998 0.086728 0.130237 0.157331 0.141638 0.114826 0.086691 0.066755 0.052660 0.042163 0.031579 0.022788 0.043606 1.000000

+ πk,1

+ πk,4

+ πk,5

+ πk,6

236

M. Yu, A. S. Alfa / European Journal of Operational Research 244 (2015) 227–239

Table 3 Joint distribution of the number of customers in the queue and the number of customers in a served batch at arbitrary epoch for DMAP/G(1,4,6)/1/12 queue. k

0¯ 0 1 2 3 4 5 6 7 8 9 10 11 12 Sum

π0¯ = π0¯ eL

r=1

r=4

r=5

r=6

Total

0.027597 0.051427 0.072371 0.044222 0.022098 0.007924 0.001539 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.227178

0.037433 0.031937 0.026889 0.021892 0.016955 0.012195 0.007862 0.004313 0.001863 0.000557 0.000085 0.000000 0.000000 0.161981

0.030424 0.026344 0.022592 0.018872 0.015174 0.011547 0.008108 0.005065 0.002673 0.001104 0.000313 0.000045 0.000000 0.142261

0.023047 0.037930 0.047105 0.051948 0.052769 0.050381 0.052956 0.042863 0.033405 0.024853 0.017562 0.011743 0.015403 0.461965

0.006615 0.118501 0.147638 0.168957 0.136934 0.106996 0.082047 0.070465 0.052241 0.037941 0.026514 0.017960 0.011788 0.015403 1.000000

 πk,1

 πk,4

 πk,5

 πk,6

0.006615

0.006615

Table 4 Joint distribution of the number of customers in the queue and the number of customers in a served batch at pre-arrival epoch for D-MAP/G(1,4,6)/1/12 queue. k

0¯ 0 1 2 3 4 5 6 7 8 9 10 11 12 Sum

π0¯−

r=1

r=4

r=5

r=6

Total

0.027160 0.050970 0.072165 0.044598 0.022384 0.008062 0.001575 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.226914

0.036996 0.031981 0.026974 0.021981 0.017041 0.012276 0.007930 0.004362 0.001891 0.000567 0.000087 0.000000 0.000000 0.162086

0.030103 0.026377 0.022655 0.018938 0.015240 0.011610 0.008165 0.005111 0.002704 0.001120 0.000319 0.000047 0.000000 0.142389

0.022794 0.037749 0.047000 0.051918 0.052801 0.050457 0.053004 0.043024 0.033563 0.024990 0.017672 0.011827 0.015550 0.462349

0.006262 0.117053 0.147077 0.168794 0.137435 0.107466 0.082405 0.070674 0.052497 0.038158 0.026677 0.018078 0.011874 0.015550 1.000000

− πk,1

− πk,4

− πk,5

− πk,6

0.006262

0.006262

When the inter-arrival time follows geometric distribution, the joint distribution of the number of customers in the queue and the number of customers in a served batch at pre-arrival and arbitrary epochs should be identical due to the memoryless property of Bernoulli arrivals. By numerical examples, we show that such a relationship is indeed established in batch-size-dependent service queue. Example 2. The numerical results are still given for the DMAP/G(1,4,6)/1/12 queue. The service time distribution for each batch size r is the same as that of Example 1. We take arrival matrices     as C = 0.3 and D = 0.7 . This means that the expected number of arrivals per unit time is also 0.7, and the discrete-time Markovian arrival process becomes the simplest point process, namely Bernoulli process, for which the matrices C and D are scalars. It can be seen from Table 5 that pre-arrival and arbitrary epoch  = π − ) due to probabilities are the same (i.e. π ¯ = π ¯− and πk,r k,r 0 0 Bernoulli arrivals. It is obvious that the numerical experiment result coincides with the BASTA (i.e. Bernoulli arrivals see time averages) property and indicates the correctness of the theoretical analysis in the previous section and the reliability of the computation program. Therefore, this method often acts as an effective internal check on our results. Example 3. Having demonstrated the applicability of the analytical results in tabular form, we have a look at how the distributions and the variances of service time affect the queueing process. In this example, numerical results are given for the DMAP/G(1,5,9)/1/N queue by varying N from 15 to 40. The input parameters of DMAP are chosen as

&

C=

0.15 0.1

0.15 0.1

'

&

,

D=

0.35 0.3

0.35 0.5

'

,

which leads to λ∗ = 0.76. We focus our attention on the loss probability of an arbitrary customer and the average waiting time of a customer in the queue. Toward this end, we first take three sets of batch-size-dependent service time as follows: Case 1.

1 1 , j = 1, 2, . . . , 5, Pr {S5 = j} = , 5 13 j = 1, 2, . . . , 13, 1 1 Pr {S6 = j} = , j = 1, 2, . . . , 15, Pr {S7 = j} = , 15 17 j = 1, 2, . . . , 17, 1 1 Pr {S8 = j} = , j = 1, 2, . . . , 19, Pr {S9 = j} = , 19 21 j = 1, 2, . . . , 21.

Pr {S1 = j} = service at arbitrary and pre-arrival epochs are displayed in the last column and last row of Tables 3 and 4, respectively. Here, all the numerical results are presented without phases. The following numerical example is useful in debugging Matlab programs and checking accuracy for computations. Once the stationary distribution π is obtained, the joint distribution of the number of customers in the queue and the number of customers in a served batch − , can be computed as well. at pre-arrival epoch, namely π ¯− and πk,r 0

Table 5 Joint distribution of the number of customers in the queue and the number of customers in a served batch at arbitrary and pre-arrival epochs for DMAP/G(1,4,6)/1/12 queue. k

0¯ 0 1 2 3 4 5 6 7 8 9 10 11 12

π0¯

0.005926

π0¯−

r=1

r=4

r=5

r=6

 πk,1

− πk,1

 πk,4

− πk,4

 πk,5

− πk,5

 πk,6

− πk,6

0.026734 0.050802 0.072550 0.044738 0.022238 0.007826 0.001452 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000

0.026734 0.050802 0.072550 0.044738 0.022238 0.007826 0.001452 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000

0.037469 0.032406 0.027344 0.022289 0.017273 0.012405 0.007939 0.004275 0.001781 0.000501 0.000070 0.000000 0.000000

0.037469 0.032406 0.027344 0.022289 0.017273 0.012405 0.007939 0.004275 0.001781 0.000501 0.000070 0.000000 0.000000

0.030409 0.026655 0.022901 0.019149 0.015410 0.011724 0.008204 0.005071 0.002615 0.001036 0.000275 0.000036 0.000000

0.030409 0.026655 0.022901 0.019149 0.015410 0.011724 0.008204 0.005071 0.002615 0.001036 0.000275 0.000036 0.000000

0.022921 0.037969 0.047339 0.052228 0.052995 0.050527 0.052682 0.042676 0.033169 0.024562 0.017252 0.011460 0.014716

0.022921 0.037969 0.047339 0.052228 0.052995 0.050527 0.052682 0.042676 0.033169 0.024562 0.017252 0.011460 0.014716

0.005926

M. Yu, A. S. Alfa / European Journal of Operational Research 244 (2015) 227–239

237

Table 6 PH representation for service time of batch size r. r

βr

Tr

Case 1 1 6 8



1  e 5 5



1  e 15 15



1  e 19 19

Case 2  3  1 e3 6 8

 

 5  e5 

 e10 10

Case 3  3  1 e3 6 8

 

 5  e5 

e10 10







E [Sr ]

01×4 0 I 4 04×1

01×9 0 I 9 09×1 01×2 0 I 2 02×1 01×7 0 I 7 07×1 01×9 0 I 9 09×1

βr

5

1  e 13 13

Tr 

3

01×18 0 I 18 018×1

01×7 0 I 7 07×1

r



01×14 0 I 14 014×1

01×2 0 I 2 02×1

Var [Sr ]

2

 8

18.67



1  e 17 17

7

 10

30

9

1  e 21 21

3

0

5

1  e 13 13







 8

0

1  e 17 17

7

 10

0

9

1  e 21 21

3

0

5

 7  e7

7

 9  e9

  8

0

 10

0

 

  

 11  e11

9

E [Sr ]

01×12 0 I 12 012×1 01×16 0 I 16 016×1 01×20 0 I 20 020×1 01×12 0 I 12 012×1 01×16 0 I 16 016×1 01×20 0 I 20 020×1 01×6 0 I 6 06×1 01×8 0 I 8 08×1

Var [Sr ]

 7

14

9

24

  11

36.67

 7

14

9

24

  11

36.67

 7

0

9

0

11

0



01×10 0 I 10 010×1



Case 2.

Pr {S1 = 3} = 1, Pr {S5 = j} = Pr {S6 = 8} = 1, Pr {S7 = j} = Pr {S8 = 10} = 1, Pr {S9 = j} =

1 , 13 1 , 17 1 , 21

j = 1, 2, . . . , 13, j = 1, 2, . . . , 17, j = 1, 2, . . . , 21.

0.05

0.04

Case 3.

Pr {S1 = 3} = 1, Pr {S5 = 7} = 1, Pr {S6 = 8} = 1, Pr {S7 = 9} = 1, Pr {S8 = 10} = 1, Pr {S9 = 11} = 1.

Ploss

0.03

0.02

Their corresponding PH representations are given in Table 6.

Example 4. In this example, we show how the performance characteristics of our model differ from general PH renewal process which has independent arrivals. For simplicity, we take a special case of DMAP, namely the 2-state Markov modulated Bernoulli process (MMBP). A 2-state MMBP is characterized by its transition probabilˆ of arrival probabilities, as given ity matrix Pˆ and a diagonal matrix  below,

Pˆ =

p 1−q

1−p q

'

and

ˆ = 

&

λ1 0

0

λ2

'

0.01

0

−0.01 15

20

25

30

35

Fig. 1. Effect of N on Ploss for different service time distributions.

16

14

Case 1 Case 2 Case 3

12

10

8

6

4 15

20

25

30

35

N

,

40

N

Wq

As we have seen in Table 6, for three different cases, the corresponding service time distributions are characterized by the same mean value. However, these are qualitatively different in that they have different variances. The total variances of above three groups of service time distribution are 125.34, 74.67 and 0, respectively. Fig. 1 displays the impact of buffer size on the loss probability for various service time distributions. We see that Case 3 has the lowest loss probability and remains almost zero when N changes from 15 to 40, while in Case 1 and Case 2, Ploss decreases exponentially as the buffer size increases. It can be readily observed that the loss probability of the system worsens as the variances of the service time distributions become larger, and the decreasing rate of the loss probability becomes slower as the buffer size increases. This corresponds with our intuitive expectations. The effect of buffer size on average waiting time of a customer in the queue is studied in Fig. 2. We can observe from this figure that for fixed N, the value of Wq is much lower in Case 3 as compared to Cases 1 and 2. Wq remains almost constant in the third case as N increases. These analytical results obtained from above numerical experiments lead to a conclusion that variability of service times does affect main performance measures even though they have the same mean value.

&

Case 1 Case 2 Case 3

Fig. 2. Effect of N on Wq for different service time distributions.

40

238

M. Yu, A. S. Alfa / European Journal of Operational Research 244 (2015) 227–239 Table 7 PH representation for service time distribution of batch size r. r

βr

1

1  e 11 11

6

1  e 17 17

Tr   01×10 0 010×1 I  10  01×16 0 I 16 016×1

E [Sr ]

Var [Sr ]

r

βr

6

10

5

1  e 15 15

9

24

7

1  e 19 19

where 0 < p, q < 1 and 0 < λ1 , λ2 < 1. The total arrival rate, λ∗ , the squared coefficient of variation, CV2MMBP(2), of inter-arrival times between any two consecutive customers, and the lag-1 autocorrelation coefficient of inter-arrival times, φMMBP(2)(1), are given respectively as follows (see, Goswami & Selvaraju, 2010):

λ1 (1 − q) + λ2 (1 − p) 2−p−q

λ1 (1 − q) + λ2 (1 − p)

2−p−q φMMBP(2)(1) =

CV2MMBP(2)

− 1,

1 1 , j = 1, 2, . . . , 11, Pr {S5 = j} = , 11 17 j = 1, 2, . . . , 17, 1 1 Pr {S6 = 8} = , j = 1, 2, . . . , 19, Pr {S7 = j} = , 19 21 j = 1, 2, . . . , 21.

Their corresponding PH representations are given in Table 7. Furthermore, the arrival process coded as MMBP(2) is defined by two matrices

ˆ = 

&

0.9814815 0.0185185

0.0185185 0.9814815

0.5000000 0.0000000

0.0000000 0.1000000

'

and ' .

Thus, the corresponding DMAP representation of above MMBP process is given by

 &  ˆ = 0.4907405 C MMBP(2) = Pˆ I 2 −  0.0092595 and

ˆ = DMMBP(2) = Pˆ 

&

0.4907405 0.0092595

0.0166671 0.8833329

0.0018519 0.0981481

'

' .

Meanwhile, the arrival process coded as PH is defined by the vector   α = 0.8333333, 0.1666667 and the matrix

&

M=

0.4907405 0.0092595

10

30

0.008

0.006 0.005 0.004 0.003 0.002

M M BP(2), λ∗ = 0.3

0.001

P H, λ∗ = 0.3

λ1 λ2 (λ1 −λ2 )2 (1−p)(1−q)(p+q−1)2 . (2−p−q)2 [λ1 (1−q)+λ2 (1−p)+λ1 λ2 (p+q−1)]2

Pr {S1 = j} =

&

18.67

0.007

The squared coefficient of variation, CV2MMBP(2), of the inter-arrival times is an important measure of the degree of traffic burstiness which is defined as the ration of the variance to the square of mean interarrival times. The lag-1 autocorrelation coefficient, which is the ratio of covariance of inter-arrival times to the variance, is used to tell us how much correlation is there in the arrival process. Next, the comparison results for MMBP(2)/G(1,5,7)/1/N and PH/G(1,5,7)/1/N queues are presented by considering the following batch-size-dependent service time distributions

Pˆ =

8

Var [Sr ]

0.009

,

   ! 2 λ1 1 − q + λ2 1 − p       CV2MMBP(2) = λ1 1 − q + λ2 1 − p + λ1 λ2 p + q − 1    !    !   2 λ1 1 − p + λ2 1 − q λ1 1 − q + λ2 1 − p p + q − 1 +  2      ! 2 − p − q λ1 1 − q + λ2 1 − p + λ1 λ2 p + q − 1 −

E [Sr ]

0.01

P loss

λ∗ =

Tr   01×14 0 014×1 I  14  01×18 0 I 18 018×1

0.0166671 0.8833329

'

.

0 12

14

16

18

20

22

24

26

N Fig. 3. Effect of N on Ploss for different correlation structures of the arrival process.

Then it can also be represented as DMAP process with

& C PH = M =

0.4907405 0.0092595

and

DPH = (e2 − Me2 )α =

&

0.0166671 0.8833329

0.4104935 0.0895063

'

0.0820989 0.0179013

' .

All these special DMAP processes are characterized by the same total arrival rate and same squared coefficient of variation equal to 0.3 and 2, respectively. However, these processes are qualitatively different in that they have different correlation structures. The arrival process labeled MMBP(2) has correlated arrivals with correlation between two successive inter-arrival times given by φMMBP(2)(1) = 0.264063. The second arrival process, namely PH, corresponds to renewal process and so the autocorrelation is 0. In Fig. 3, the difference between loss probabilities of these MMBP(2) and PH models is shown for various values of buffer capacity N. The numerical results indicate that the loss probability of a batch-size-dependent service model with a correlation structure cannot be assumed to be identical with that having independent arrivals, especially when the capacity of the buffer space is not big enough. Additionally, the effect of buffer size on average sojourn time of a customer in the system is examined in Fig. 4. It can be seen from this figure that Ws has higher value in MMBP(2)/G(1,5,7)/1/N model as compared to PH/G(1,5,7)/1/N model. We also reveal that the mean sojourn time in the PH/G(1,5,7)/1/N queue is not significantly affected by the changes of N. From these observations, we can conclude that there are some obvious differences in loss probabilities and mean sojourn times when we consider the correlation in the arrival process to that of the case where the arrivals are assumed to be independent. The autocorrelated arrivals often have a profound impact on system operation and usually cause the worse of system performance. Therefore, it frequently leads to inaccurate predication of performance measures under the assumption of independent arrivals.

M. Yu, A. S. Alfa / European Journal of Operational Research 244 (2015) 227–239

16.8

16.6

Ws

16.4

16.2

16 M M BP(2), λ∗ = 0.3 15.8

15.6 12

P H, λ∗ = 0.3

14

16

18

20

22

24

26

N Fig. 4. Effect of N on Wq for different correlation structures of the arrival process.

5. Conclusions In this paper, we considered a discrete-time batch-sizedependent service model and associated algorithmic aspects for a DMAP/G(1,a,b)/1/N system. Notice that any distribution with finite support on the nonnegative integers can be expressed as a discrete PH distribution, an efficient and simple procedure for the computation of the joint state probabilities at various time epochs is developed. Applying the UL-type RG factorization, the numerical experiments demonstrate that the proposed algorithm well suits for computer realization even for large number of system states. Hence, the method used in this paper seems to be applicable to a much wider range of discrete-time finite-buffer queueing system with more general batch arrival and bulk service mechanisms. Here, we need to point out that using the approach presented in the paper of Banerjee and  and π − , Gupta (2012), we cannot obtain the joint probabilities πN,r N,r r = 1, a, . . . , b. Thus, some key performance indicators such as the loss probability and the mean sojourn time of a customer are unable to derive. As a further research, it will be more interesting to analyze the actual waiting time distribution of an arbitrary customer in the system and then obtain its moments of higher orders. Acknowledgments The authors would like to thank the anonymous reviewers for their valuable comments and suggestions to improve the quality of this paper. This research was partially supported by grant from NSERC DAS programs, National Natural Science Foundation of China, China (Nos. 71301111, 71171138, 71402072) and the FSUSE (No. 2012RC23). References Alfa, A. S. (2004). Markov chain representations of discrete distributions applied to queueing models. Computers & Operations Research, 31(14), 2365–2385.

239

Alfa, A. S., & He, Q. M. (2008). Algorithmic analysis of the discrete time GIX /GY /1 queueing system. Performance Evaluation, 65(9), 623–640. Alfa, A. S., & Li, W. (2001). Matrix-geometric analysis of the discrete time GI/G/1 system. Stochastic Models, 17(4), 541–554. Banerjee, A., & Gupta, U. C. (2012). Reducing congestion in bulk-service finitebuffer queueing system using batch-size-dependent service. Performance Evaluation, 69(1), 53–70. Banerjee, A., Gupta, U. C., & Goswami, V. (2014). Analysis of finite-buffer discrete-time batch-service queue with batch-size-dependent service. Computers & Industrial Engineering, 75, 121–128. Banik, A. D., Chaudhry, M. L., & Gupta, U. C. (2008). On the finite buffer queue with renewal input and batch markovian service process: GI/BMSP/1/N. Methodology and Computing in Applied Probability, 10(4), 559–575. Bar-Lev, S. K., Parlar, M., Perry, D., Stadje, W., & Schouten, F. A. V. (2007). Applications of bulk queues to group testing models with incomplete identification. European Journal of Operational Research, 183(1), 226–237. Chakravarthy, S. R. (1992). A finite capacity GI/PH/1 queue with group services. Naval Research Logistics, 39(3), 345–357. Chakravarthy, S. R. (1993). Analysis of a finite MAP/G/1 queue with group services. Queuing Systems, 13(4), 385–407. Chang, S. H., & Choi, D. W. (2005). Performance analysis of a finite-buffer discretetime queue with bulk-arrival, bulk-service and vacations. Computers & Operations Research, 32(9), 2213–2234. Chaudhry, M. L., & Chang, S. H. (2004). Analysis of the discrete-time bulk-service queue Geo/GY /1/N+B. Operation Research Letters, 32(4), 355–363. Chaudhry, M. L., & Gai, J. (2012). A simple and extended computational analysis of M/G(j a,b)/1 and M/G(j a,b)/1/(B+b) queues using roots. INFOR, 50(2), 72–79. Chaudhry, M. L., & Gupta, U. C. (2003). Analysis of a finite-buffer bulk-service queue with discrete-markovian arrival process: D-MAP/Ga,b /1/N. Naval Research Logistics, 50(4), 345–363. Claeys, D., Steyaert, B., Walraevens, J., Laevens, K., & Bruneel, H. (2013). Analysis of a versatile batch-service queueing model with correlation in the arrival process. Performance Evaluation, 70(4), 300–316. Claeys, D., Walraevens, J., Laevens, K., & Bruneel, H. (2010). A queueing model for general group screening policies and dynamic item arrivals. European Journal of Operational Research, 207(2), 827–835. Curry, G. L., & Feldman, R. M. (1985). An M/M/1 queue with a general bulk service rule. Naval Research Logistics Quarterly, 32(4), 595–603. Goswami, C., & Selvaraju, N. (2010). The discrete-time MAP/PH/1 queue with multiple working vacations. Applied Mathematical Modelling, 34(4), 931–946. Gupta, U. C., & Goswami, V. (2002). Performance analysis of finite buffer discretetime queue with bulk service. Computers & Operations Research, 29(10), 1331–1341. Gupta, U. C., & Laxmi, P. V. (2001). Analysis of the MAP/Ga,b /1/N queue. Queueing System, 38(2), 109–124. Gupta, U. C., & Sikdar, K. (2004). The finite-buffer M/G/1 queue with general bulk-service rule and single vacation. Performance Evaluation, 57(2), 199–219. Hunter, J. (1983). Mathematical techniques of applied probability, discrete time models: Techniques and applications (Vol. 2). New York: Academic Press. Laxmi, P. V., & Gupta, U. C. (1999). On the finite-buffer bulk-service queue with general independent arrivals: GI/M[b] /1/N. Operations Research Letters, 25(5), 241–245. Li, Q. L. (2010). Constructive computation in stochastic models with applications: The RG-factorizations. New York: Springer. Li, Q. L., & Cao, J. H. (2004). Two types of rg-factorizations of quasi-birth-and-death processes and their applications to stochastic integral functionals. Stochastic Models, 20(3), 299–340. Neuts, M. F. (1967). A general class of bulk queues with poisson input. Annals of Mathematical Statistics, 38(3), 759–770. Neuts, M. F. (1975). Probability distribution of phase type. In Liber amicorum Prof. emeritus H. florin. Belgium: Department of Mathematics, University of Louvain (pp. 173–206). Samanta, S. K., Chaudhry, M. L., & Gupta, U. C. (2007). Discrete-time GeoX /G(a,b)/1/N queues with single and multiple vacations. Mathematical and Computer Modelling, 45(1–2), 93–108. Samanta, S. K., Gupta, U. C., & Sharma, R. K. (2007). Analyzing discrete-time dbmap/g/1/n queue with single and multiple vacations. European Journal of Operational Research, 182(1), 321–339. Skidar, K., & Gupta, U. C. (2005). Analytic and numerical aspects of batch service queues with single vacation. Computers & Operations Research, 32(4), 943–966. Takagi, H. (1993). Queueing analysis: A foundation of performance evaluation (Vol. 3). Amsterdam: North-Holland.