Analysis of a finite-buffer bulk-service queue under Markovian arrival process with batch-size-dependent service

Author's Accepted Manuscript

Analysis of a finite-buffer bulk-service queue under Markovian arrival process with batchsize-dependent service A. Banerjee, U.C. Gupta, S.R. Chakravarthy

www.elsevier.com/locate/caor

PII: DOI: Reference:

S0305-0548(15)00043-X http://dx.doi.org/10.1016/j.cor.2015.02.012 CAOR3733

To appear in:

Computers & Operations Research

Cite this article as: A. Banerjee, U.C. Gupta, S.R. Chakravarthy, Analysis of a finite-buffer bulk-service queue under Markovian arrival process with batchsize-dependent service, Computers & Operations Research, http://dx.doi.org/ 10.1016/j.cor.2015.02.012 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Analysis of a finite-buffer bulk-service queue under Markovian arrival process with batch-size-dependent service A. Banerjeea , U. C. Guptab,1,∗, S. R. Chakravarthyc a

Department of Mathematical Sciences, Indian Institute of Technology (BHU), Varanasi-221005. India. b Department of Mathematics, Indian Institute of Technology, Kharagpur-721302. India. c Department of Industrial and Manufacturing Engineering, Kettering University, Flint, MI-48504, USA.

Abstract We consider a finite capacity single server queue in which the customers arrive according to a Markovian arrival process. The customers are served in batches (of varying size) following a ‘general bulk service rule’. The service times, which depend on the size of the batch, are generally distributed. We obtain, in steady-state, the joint distribution of the random variables of interest at various epochs. Efficient computational procedures in the case of phase type services are presented. An illustrative numerical example to bring out the qualitative nature of the model is presented. Keywords: Algorithmic probability, General bulk service rule, Markovian arrival process, Phase type distribution, Queue 1. Introduction Bulk service queues have applications in many areas notably in production and manufacturing, transportation, package delivery, tourism, and amusement parks such as Disney World and Six Flags. For example, in production and manufacturing applications, it may not always be efficient to start the production as and when the orders arrive. Instead, one has to wait until a certain number of orders is placed to start the process. Similarly, disposing hazardous petrochemicals and petroleum wastes (that arrive in containers or drums) may need thermal treatment using high temperatures, and hence processing them efficiently require grouping them with a minimum number and maximum number of containers. In the amusement park set up people queue up for going through various thrill rides and all of these are accommodated in groups of varying sizes with restrictions on the minimum and maximum numbers in each ride. In package delivery example, it is always better to fill the trucks (of varying sizes) to their capacity to ∗

Corresponding author Email addresses: [email protected], [email protected] (A. Banerjee), [email protected] (U. C. Gupta), [email protected] (S. R. Chakravarthy) 1 Tel: +91 - 3222 - 283654

Preprint submitted to Elsevier

February 25, 2015

balance the cost/efficiency to the delivery times of the packages. Thus, one can see the variety of scenarios that call for bulk service queues in real-life applications [1–6]. Queueing systems in which the services are offered in batches of varying sizes have been extensively studied in the literature [2–18]. In such queueing systems the customers are served in batches of sizes varying from a minimum size of ‘a’ to a maximum size of ‘b’ and this service rule is referred to as the ‘general bulk service’ (GBS) rule by Neuts [16]. The book by Chaudhry and Templeton [12] provides an in-depth study of bulk service queueing systems. Thus, one can categorize the study of bulk service queueing systems as follows. The bulk service queueing systems in which the 1. buffer size is (a) finite or (b) infinite 2. arrivals occur according to a (a) renewal or (b) correlated process 3. arrivals occur (a) singly or (b) in batches 4. services are (a) independent of the batch size or (b) dependent on the size of the batch being served 5. service times are (a) exponential or (b) non-exponential. While the literature is abundant with queueing systems dealing with cases falling in the intersection of (1), (2a), (3), (4a), and (5a), very few papers deal with case (2b) in combination with 4(b) and 5(b). Earlier work on bulk service queueing systems in which the services depend on the batch size includes Curry and Feldman [19] and Neuts [20]. Further, Bar-Lev et al. [3] and Chaudhry and (a,b)

Gai [21] both considered an M/Gr

/1 queue and obtained the queue length distribution at

departure epochs. One may note here that neither Bar-Lev et al. [3] nor Chaudhry and Gai [21] obtained the queue length distribution at an arbitrary epoch which is required to compute various performance measures. Also from their analysis one cannot obtain the joint distribution of the number of customers in the queue and the number in the departing batch at a departure epoch or the joint distribution of the number of customers in the queue and the number in the service at an arbitrary epoch. These distributions are needed in order to apply any batchsize-dependent service policy to a bulk service queueing system. In view of this Banerjee and (a,b)

Gupta [22] analyzed an M/Gr

/1/N queue and obtained all such distributions. Claeys et

al. [4] analyzed an infinite buffer discrete-time bulk arrival and bulk service queue with GBS rule assuming the service time to be dependent on the size of the batch. They obtained the steady-state probability of the number in the system employing probability generating function approach. For recent development in this direction, see e.g., [23, 24], and the references therein. While the papers referenced above all dealt with renewal arrivals, there exists some literature on bulk service queueing systems [25–32] in which the customers arrive according to a versatile point process, namely, Markovian arrival process (M AP ). The M AP is a rich class of point processes that includes many well-known processes such as Poisson process, phase type (PH) renewal processes, and Markov-modulated Poisson process. It has been observed that modern 2

telecommunication networks/ ATM networks support diverse traffics with different service characteristics such as voice, data and video. In such network, IP packets or cells of voice, video, data are sent over a common transmission channel on statistical multiplexing basis. These traffic streams are statistically multiplexed and transmitted in super high speed and also they are highly irregular (e.g., bursty and correlated). A good representation of such traffic is a Markovian arrival process. Further, M AP is very useful to model situations where the inter-arrival times are no longer independent. For example, in applications where the arrivals are generated from different sources the pooled input will not necessarily follow a renewal process (even if the individual sources are renewal). This is very common in production and manufacturing, computer communications, transportation, and other areas wherein the model under study is very much applicable. When the arrivals occur in batches in the context of M AP , we refer to that arrival process as BM AP . For more details on BM AP processes and their usefulness in stochastic modelling, we refer to [33, 34], and for a review on M AP and BM AP , we refer the reader to [33, 35–39]. Chakravarthy [27] and Gupta and Laxmi [28] both analyzed single server finite buffer bulk service queues with M AP arrivals and with GBS rule for the services. Chakravarthy [27] analyzed a finite-capacity M AP/G(a,b) /1 queue where the upper threshold, ‘b’, is taken as the buffer size. Later, Gupta and Laxmi [28] considered a more general model than the one considered by Chakravarthy [27], and assumed that the buffer size, N , is greater than the upper threshold value ‘b’, and obtained the queue length distribution at various epochs. In the context of a multi-server queueing system with M AP arrivals, exponential group services with batches of sizes between two pre-determined thresholds, and a dynamic service rule that governs the size of the batch at the time of starting a service, Chakravarthy and Alfa [40] studied a bulk service queueing system using matrix-analytic methods. Recently, Banik [26] analyzed BM AP/G(a,b) /1/N and BM AP/M SP (a,b) /1/N queues and obtained the queue length distributions for both the models. But in all of these papers the service times are assumed to be independent of the serving batch size. In Chakravarthy and Alfa [41] a finite capacity queuing system with M AP arrivals attended by two exponential servers, who offer services in groups of varying sizes, is studied. Here the authors assume that the service times may depend on the number of customers in service. Efficient algorithmic procedures for computing the steady-state queue length densities and other system performance measures are discussed. A conjecture on the nature of the mean waiting time is proposed. Some illustrative numerical examples are presented. Alfa et al. [42] studied a discrete queueing system with M AP arrivals and PH services in batches of size ranging from 1 to a fixed threshold. Under a probabilistic service rule, which depends on the number of customers waiting in the queue, the authors show that the steady-state probability vector is (modified) matrix-geometric type ([43]). Using a dynamic probabilistic rule associated with group (the size of which varies from a pre-determined threshold to the maximum buffer size)

3

services, assuming M AP arrivals and exponential services whose parameter depends on the size of the batch, Chakravarthy and Bin [44] developed efficient algorithms for computing various performance measures such as throughput, mean number served, job overflow probability and server idle probability, useful in qualitative and quantitative interpretations of the model. Most of the papers dealing with batch-size-dependent services assume the service times to be exponential except in Alfa et al. [42] where it is of PH type. But this one, as mentioned earlier, is studied in discrete time. In view of these, to fill the gap in the present literature, in this paper we study a continuous time finite capacity single server queue where arrivals occur according to a M AP with representation (C, D) of order m. The services are offered in batches of size varying from a minimum size a to a maximum value of b. The service time distribution of the batch is assumed to be arbitrarily distributed and dependent on the batch size. Henceforth, we (a,b)

denote this model by M AP/Gr

/1/N . For analytic analysis of this model we use the embedded

Markov chain technique and first obtain the joint distribution of the number of customers in the queue, the number in the departing batch, and the phase of the arrival process at a departure epoch. Then using the supplementary variable technique and considering the supplementary variable as remaining service time of the batch undergoing service, we develop relations between arbitrary and departure epoch probabilities in order to obtain the joint distribution of the number of customers in the queue, the number in service, and phase of the arrival process at an arbitrary epoch. Then we obtain arrival epoch probabilities by deriving relations between arbitrary and arrival epoch probabilities. Several performance measures of interest, such as, mean number of customers in the queue (system), mean number of customers served in a batch, loss probability and mean waiting time in the queue (system), have been obtained. For use in sequel, let e(r), ej (r) and Ir denote, respectively, the column vector of dimension r consisting of 1’s, column vector of dimension r with 1 in the j th position and 0 elsewhere, and an identity matrix of dimension r. When there is no need to emphasize the dimension of these vectors we will suppress the suffix. Thus, e will denote a column vector of 1’s of appropriate dimension. The notation ‘T ’ appearing in the superscript will stand for the transpose of a matrix; the notation ‘’ will stand for the derivative; and the symbols ⊗ and ⊕, respectively, will stand for the Kronecker product and Kronecker sum of two matrices. Thus, if A is a matrix of order m × n and if B is a matrix of order p × q, then A ⊗ B will denote a matrix of order mp × nq whose (i, j)th block matrix is given by aij B; the Kronecker sum of two square matrices, say, G of order g and H of order h, is given by G ⊗ I + I ⊗ H, a square matrix of dimension gh. For more details on Kronecker products and sum, we refer the reader to Marcus and Minc [45]. The paper is organized as follows. In Section 2, the description of the model and its analysis at various epochs is given. Some key system performance measures to bring out the qualitative nature of the model are displayed in Section 3. The computational simplifications that arise when the services are of phase type are spelled out in Section 4, and an illustrative numerical example is discussed in Section 5. Some concluding remarks are presented in Section 6.

4

2. Model description and analysis We consider a single server queueing system in which the customers arrive according to a M AP where the arrivals are governed by an underlying m-state Markov chain with transition rate cij , 1 ≤ i, j ≤ m, i = j, there is a transition from state i to state j in the underlying Markov chain without an arrival, and with transition rate dij , 1 ≤ i, j ≤ m, there is a transition from state i to state j in the underlying Markov chain with an arrival. The matrix C = [cij ] has nonnegative off-diagonal and negative diagonal elements, and the matrix D = [dij ] has ˜ (t) denotes the number of customers arriving in (0, t] and J(t) be nonnegative elements. Let N the state of the underlying Markov chain at time t with state space {i : 1 ≤ i ≤ m}. Then ˜ (t), J(t)} is a two-dimensional Markov process with state space {(n, i) : n ≥ 0, 1 ≤ i ≤ m} {N and is known as Markovian arrival process (M AP ). For this arrival process we have (C+D)e = 0. As (C + D) is the infinitesimal generator of the underlying Markov chain {J(t)}, there exists a stationary probability vector δ such that δ(C + D) = 0, δe = 1.

(1)

The fundamental arrival rate (average arrival rate) of the above Markov process is given by λ∗ = δDe. For more detail on this topic see [46]. Let {P˜ (n, t), n ≥ 0, t ≥ 0} be the square matrices of order m whose (i, j)-th elements are the conditional probabilities defined as ˜ (t) = n, J(t) = j|N ˜ (0) = 0, J(0) = i}, n ≥ 0, 1 ≤ i, j ≤ m. p˜i,j (n, t) = P {N ˜ (t), J(t); t ≥ 0}, satisfy the following These matrices, associated with the counting process {N system of difference-differential equations P˜  (0, t) = P˜ (0, t)C, P˜  (n, t) = P˜ (n, t)C + P˜ (n − 1, t)D, with

n ≥ 1,

P˜ (0, 0) = I.

and have been extensively studied in the literature with an efficient procedure for computing them is given in [47]. There is a finite waiting room of size N such that any arrival finding the buffer full will be considered lost. Suppose that a and b satisfying 1 ≤ a ≤ b < N are two pre-determined thresholds. The server offers services in batches of varying sizes according to the GBS rule. That is, if the number of customers waiting in the queue at a service completion is less than a, then the server waits until the number waiting hits a and then offers services to all a customers in a batch; otherwise, the server will take a batch of r, a ≤ r ≤ b, customers using first-come-firstserved rule and offers services to the entire batch. The service times are assumed to be generally distributed and may depend on the batch size. Specifically, let Tr , a ≤ r ≤ b, denote the service 5

˜ r (.) and H ∗ (.), a ≤ r ≤ b, time for a batch of size r with distribution function Hr (.). Let h r denote, respectively, the probability density function and the Laplace-Stieltjes transform of Tr . Let hr denote the mean service time for a batch of size r. The steady-state analysis of the model under study will be carried out using the embedded Markov chain approach since the service times are assumed to be generally distributed. First, we will look at the semi-Markov process embedded at points of departure of (batches) of customers. Towards this end, we define the following conditional probabilities. (r)

• [Ak (x)]ij , a ≤ r ≤ b, k ≥ 0, 1 ≤ i, j ≤ m, x ≥ 0, is the conditional probability that, starting with a departure which left at least a customers in the queue with the arrival process in state i, the next departure occurs no later than time x, and during the service period of r customers exactly k new customers arrive; the phase of the arrival process is in phase j at the departure epoch. • [Bn,k (x)]ij , 0 ≤ n ≤ a − 1, k ≥ 0, 1 ≤ i, j ≤ m, x ≥ 0, is the conditional probability that, starting with a departure which left n customers in the queue with the arrival process in state i, the next departure occurs no later than time x, and and during the service period of a customers exactly k new customers arrive; the phase of the arrival process is in phase j at the departure epoch. (r)

Denote by Ak (x), a ≤ r ≤ b, k ≥ 0, x ≥ 0, and Bn,k (x), 0 ≤ n ≤ a − 1, k ≥ 0, x ≥ 0, respectively, (r)

the square matrices of order m, whose (i, j)th elements are given by [Ak (x)]ij and [Bn,k (x)]ij . (r)

Hence using the definition of Ak (x), Bn,k (x), we have (r)

Ak (x) =

x 0

P˜ (k, t)dHr (t), 0 ≤ k ≤ N, a ≤ r ≤ b,

(b) A˜k (x) 

(r) A˜N +1 (x) 

Bn,k (x) =

x 0

∞  j=k

(b)

Aj (x), b ≤ k ≤ N,

∞  j=N +1

(r)

Aj (x), a ≤ r ≤ b − 1,

(a) P˜ (a − 1 − n, x − u)DAk (u)du, 0 ≤ n ≤ a − 1, 0 ≤ k ≤ N,

˜n,N +1 (x)  B

∞  k=N +1

Bn,k (x), 0 ≤ n ≤ a − 1 .

6

For use in sequel, we define ˜  (−C)−1 D, D (r)

(r)

Ak  Ak (∞), 0 ≤ k ≤ N, a ≤ r ≤ b, (b) (b) (r) (r) A˜k  A˜k (∞), b ≤ k ≤ N, A˜N +1  A˜N +1 (∞), a ≤ r ≤ b − 1,

(2)

˜ a−n A(a) , 0 ≤ n ≤ a − 1, 0 ≤ k ≤ N, Bn,k  Bn,k (∞) = D k ˜n,N +1 (∞) = D ˜ a−n A˜(a) , 0 ≤ n ≤ a − 1. ˜n,N +1  B B N +1 ∞ ˜ n+1 , n ≥ 0, [28]. The last two equations in (2) follow from the fact that 0 P˜ (n, t)Ddt = D 2.1. Joint distribution at departure epoch In this section, we obtain the the joint distribution of the number of customers in the queue, the number in the departing batch, and the phase of the arrival process at a departure epoch. Towards this end, we first define, for n ≥ 1, the following variables. • Nn+ denotes the number of customers in the queue at the nth departure epoch. • Sn+ denotes the number of customers in the departing batch at the epoch of the nth departure. • Jn+ denotes the phase of the arrival process immediately after the nth departure. The process {(Nn+ , Sn+ , Jn+ ) : n ≥ 1}, constitutes a three-dimensional Markov chain in discretetime with state space {(i, r, k): 0 ≤ i ≤ N, a ≤ r ≤ b, 1 ≤ k ≤ m} with transition probability matrix (P ) of dimension (N + 1)(b − a + 1)m given by

0

⎛

0

1

...

N −b

...

N −1

B0,0

B0,1

...

B0,N −b

...

B0,N −1

B1,1 .. .

... .. .

B1,N −b .. .

... ...

B1,N −1 .. .

Ba−1,1

...

Ba−1,N −b

...

Ba−1,N −1

...

AN −1

...

AN −1 .. .

⎜ ⎜ B1,0 ⎜ ⎜ .. ⎜ . ⎜ ⎜ a − 1 ⎜ Ba−1,0 ⎜ ⎜ A(1) a ⎜ 0 ⎜ (2) ⎜ a + 1 A P = 0 ⎜ ⎜ .. .. ⎜ . . ⎜ ⎜ (b−a+1) b ⎜ A0 ⎜ b+1⎜ 0 ⎜ ⎜ .. .. ⎜ . . ⎝ N 0 1 .. .

(1) A1 (2) A1

.. .

(b−a+1) A1 (b−a+1) A0

... ... .. . ...

.. .

... .. .

0

...

(1) AN −b (2) AN −b

.. .

(b−a+1) AN −b (b−a+1) AN −b−1

...

(b−a+1) A0

7

...

(1) (2)

(b−a+1)

...

AN −1

... .. .

AN −2 .. .

...

Ab−1

(b−a+1)

(b−a+1)

N B¯0,N

⎞

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ B¯a−1,N ⎟ ⎟ (1) ⎟ A¯N ⎟ ⎟ (2) ⎟ A¯N ⎟ ⎟ .. ⎟ . ⎟ ⎟ (b−a+1) A¯N ⎟ (b−a+1) ⎟ ¯ AN −1 ⎟ ⎟ ⎟ ⎟ ⎠ (b−a+1) A¯ B¯1,N .. .

b

Each ‘0’ in the above matrix is a null matrix of order (b − a + 1)m × (b − a + 1)m. The other (l) block elements of P are Bi,j ’s (0 ≤ i ≤ a − 1, 0 ≤ j ≤ N − 1), B¯i,N ’s (0 ≤ i ≤ a − 1), A ’s j

(b−a+1) ’s (b ≤ j ≤ N ) are matrices of dimension (0 ≤ j ≤ N − 1, 1 ≤ l ≤ b − a + 1) and A¯j

(b − a + 1)m × (b − a + 1)m and are given as follows:

1

Bi,j

⎛

1

(l)

Aj

Bi,N

0 ...

0

...

0 ... .. . . . .

0 .. .

... .. .

0 ... .. . . . .

0 .. .

... .. .

0 ...

0

...

...

ζ

Bi,j

0

...

0

...

0

0 .. .

... .. .

0 .. .

... .. .

0 ... .. . . . .

0 .. .

... .. .

0

0

...

⎟ ⎜ 0⎟ 2⎜ ⎟ ⎜ Bi,N .. ⎟ .. ⎜ .. ⎜ .⎟ ⎟, B¯i,N = . ⎜ . ⎟ 0⎟ l⎜ ⎜ Bi,N ⎟ .. ⎟ .. ⎜ .. .⎠ .⎜ ⎝ . 0 ζ Bi,N

...

ζ

...

0

...

1 2 ... 0 0 ...

0 0 ...

l (r) Aj (r) Aj

⎞

1

⎛

1 2

...

0 0

. . . AN

...

. . . AN .. .. . .

... .. .

. . . AN .. .. . .

... .. .

. . . AN

...

.. .

Aj .. .

... .. .

Aj

...

0

(r)

(r)

...

0 0 ... ⎜ 2⎜ 0 0 . . . ⎜ .. ⎜ .. .. . . . .⎜ . . = ⎜ ⎜ l ⎜0 0 ... .. ⎜ ⎜ . . .. . . ⎝ .. ..

0

...

0 .. .

... .. .

0 .. .

... .. .

0

...

ζ

1

⎛

⎟ ⎜ 0⎟ 2⎜ 0 0 ⎟ ⎜ .. ⎜ .. .. .. ⎟ ⎜ . .⎟ ⎟, A(l) = ⎜ . . N ⎟ 0⎟ l⎜ ⎜0 0 .. ⎜ .. ⎟ ⎟ ⎜. . . ⎝ .. .. .⎠

l

⎛

⎞

... .. .

1 2 ...

1

(b−a+1) A¯j

...

l

⎜ 2⎜ 0 0 . . . ⎜ .. ⎜ .. .. . . . .⎜ . . = ⎜ ⎜ l ⎜0 0 ... .. ⎜ ⎜ . . .. . . ⎝ .. .. ζ

l

...

Bi,j

⎛

2 ...

2

⎜ 2⎜ ⎜ Bi,j .. ⎜ .. .⎜ . = ⎜ l⎜ ⎜ Bi,j .. ⎜ .. .⎜ ⎝ . ζ

1

1

0 0 ...

ζ ζ (b) ⎞ A˜j (b) ⎟ A˜j ⎟ ⎟ .. ⎟ . ⎟ ⎟ (b) A˜j ⎟ ⎟ .. ⎟ ⎟ . ⎠ (b) A˜

0 0

l (r) (r)

(r)

(r)

...

ζ

⎞ ˜i,N +1 B ⎟ ˜i,N +1 ⎟ B ⎟ ⎟ .. ⎟ . ⎟, ˜ Bi,N +1 ⎟ ⎟ ⎟ .. ⎟ . ⎠ ˜i,N +1 B ζ

(r) ⎞ A˜N +1 (r) ⎟ A˜N +1 ⎟ ⎟ .. ⎟ . ⎟ ⎟, (r) A˜N +1 ⎟ ⎟ .. ⎟ ⎟ . ⎠ (r) A˜N +1

for b ≤ j ≤ N ,

j

where r = l + a − 1, each elememt ‘0’ of the above matrices is a m × m zero matrix and ζ = b − a + 1 is a constant. Let π + = (π + (0), · · · , π + (N )) denote the steady-state probability vector of P . That is, π + P = π + , π + e = 1.

(3)

+ th component of We further partition π + (n) = (π + a (n), · · · , π b (n)), 0 ≤ n ≤ N . Note that the j

the vector, π + r (n), of dimension m, gives the steady-state probability that at a departure epoch there are n customers in the system, the departing batch size is r, and at that time the phase of the arrival process is in j, for 0 ≤ n ≤ N, a ≤ r ≤ b, 1 ≤ j ≤ m. 8

2.2. Joint distribution at an arbitrary epoch In this section, we obtain the joint distribution of the number of customers in the queue, the number in service, and phase of the arrival process at an arbitrary epoch. Towards this end, we define the following quantities at an arbitrary epoch t. • Nq (t) is the number of customers in the queue waiting for service • S(t) is the number of customers in service, if any • J(t) is the phase of the arrival process • U (t) is the remaining service time (of a batch), if any Let us define the following joint probabilities, for i, 1 ≤ i ≤ m, pˆi (n, 0; t) = P {Nq (t) = n, S(t) = 0, J(t) = i}, 0 ≤ n ≤ a − 1, π ˆi (n, r, u; t)du = P {Nq (t) = n, S(t) = r, J(t) = i, u < U (t) ≤ u + du}, 0 ≤ n ≤ N, a ≤ r ≤ b, u ≥ 0, and in limiting case 1 ≤ i ≤ m, 0 ≤ n ≤ a − 1,

pi (n, 0) = lim pˆi (n, 0; t), t→∞

ˆi (n, r, u; t), π ˜i (n, r, u) = lim π t→∞

1 ≤ i ≤ m, 0 ≤ n ≤ N, a ≤ r ≤ b.

˜ (n, r, u) = (˜ π1 (n, r, u), · · · , π ˜m (n, r, u)). Let us define the vectors p(n, 0) = (p1 (n, 0), · · · , pm (n, 0)) and π Now relating the state of the system at time t and t+dt, by considering each phase and using the supplementary variable technique, in steady-state, we have the following differential equations in matrix form. 0 = p(0, 0)C +

b

π ˜ (0, j, 0),

(4)

π ˜ (n, j, 0), 1 ≤ n ≤ a − 1,

(5)

j=a

0 = p(n, 0)C + p(n − 1, 0)D +

b

j=a

b

−

d ˜ a (u) + ˜ a (u), π ˜ (0, a, u) = π ˜ (0, a, u)C + p(a − 1, 0)D h π ˜ (a, j, 0)h du

(6)

j=a

b

d ˜ r (u), a + 1 ≤ r ≤ b , ˜ (0, r, u) = π ˜ (0, r, u)C + π ˜ (r, j, 0)h − π du

(7)

j=a

−

d π ˜ (n, r, u) = π ˜ (n, r, u)C + π ˜ (n − 1, r, u)D, 1 ≤ n ≤ N − 1, a ≤ r ≤ b − 1, du

9

(8)

b

−

d ˜ b (u), 1 ≤ n ≤ N − b, π ˜ (n + b, j, 0)h π ˜ (n, b, u) = π ˜ (n, b, u)C + π ˜ (n − 1, b, u)D + du

(9)

j=a

d π ˜ (n, b, u) = π ˜ (n, b, u)C + π ˜ (n − 1, b, u)D, N − b + 1 ≤ n ≤ N − 1, du d ˜ (N, r, u) = π ˜ (N, r, u)(C + D) + π ˜ (N − 1, r, u)D, a ≤ r ≤ b. − π du For detail derivation of (4) - (11) see Appendix A. −

(10) (11)

˜ r (u) as Define the Laplace transforms of π ˜ (n, r, u) and h ∞ ˜ (n, r, u)du, a ≤ r ≤ b, 0 ≤ n ≤ N, Re s ≥ 0, π ∗ (n, r, s) = 0 e−su π Hr∗ (s) =

∞ 0

˜ r (u)du, a ≤ r ≤ b, Re s ≥ 0, e−su h

and observe that ∞

∗

π ˜ (n, r, u)du, a ≤ r ≤ b, 0 ≤ n ≤ N.

π(n, r)  π (n, r, 0) =

(12)

0

Note that the ith component of the vector, π(n, r), of dimension m, gives the steady-state probability that at arbitrary-epoch there are n customers in the queue, the server is busy in serving a batch of size r, and at that time the phase of the arrival process is in i, for 0 ≤ n ≤ N, a ≤ r ≤ b, 1 ≤ i ≤ m. Now multiplying equations (6-11) by e−su and integrating with respect to u over 0 to ∞, we have ˜ (0, a, 0) = π ∗ (0, a, s)C + Ha∗ (s) −sπ ∗ (0, a, s) + π

b

π ˜ (a, j, 0) + Ha∗ (s)p(a − 1, 0)D,

(13)

j=a

˜ (0, r, 0) = π ∗ (0, r, s)C + Hr∗ (s) −sπ ∗ (0, r, s) + π

b

π ˜ (r, j, 0), a + 1 ≤ r ≤ b,

(14)

j=a

˜ (n, r, 0) = π ∗ (n, r, s)C + π ∗ (n − 1, r, s)D, a ≤ r ≤ b − 1, 1 ≤ n ≤ N − 1, (15) −sπ ∗ (n, r, s) + π ˜ (n, b, 0) = π ∗ (n, b, s)C +Hb∗ (s) −sπ ∗ (n, b, s)+ π ∗

∗

b

π ˜ (n+b, j, 0)+π ∗ (n−1, b, s)D, 1 ≤ n ≤ N −b,

j=a

(16)

∗

˜ (n, b, 0) = π (n, b, s)C + π (n − 1, b, s)D, N − b + 1 ≤ n ≤ N − 1, −sπ (n, b, s) + π

(17)

˜ (N, r, 0) = π ∗ (N, r, s)(C + D) + π ∗ (N − 1, r, s)D, a ≤ r ≤ b. −sπ ∗ (N, r, s) + π

(18)

Post-multiplying (4) and (5) by e and adding the resulting equations, and using the fact that (C + D)e = 0, we get p(a − 1, 0)De =

b a−1

n=0 r=a

10

π ˜ (n, r, 0)e.

(19)

Similarly, post-multiplying (13-18) by e and adding the resulting equations, and noting that (C + D)e = 0, we get b N

n=0

b b a b



1 − Ha∗ (s)

1 − Hn∗ (s) π ˜ (n, r, 0)e π (n, r, s)e = π ˜ (n, r, 0)e + s s r=a r=a r=a ∗

n=0

n=a+1

(20) +

1 − Hb∗ (s) s

N

b

π ˜ (n, r, 0)e.

n=b+1 r=a

Letting s → 0 in (20) and using the fact that a−1

p(n, 0) +

b N



π(n, r) = δ,

(21)

n=0 r=a

n=0

where δ is as given in (1), we obtain 1−

a−1

p(n, 0)e = ha

b a



π ˜ (n, r, 0)e +

n=0 r=a

n=0

b b



N b



hn π ˜ (n, r, 0)e + hb

n=a+1 r=a

π ˜ (n, r, 0)e. (22)

n=b+1 r=a

Before we establish the result relating the probability vectors {p(n, 0), π(n, r)} and {π + (n, r), ˜ (n, r, 0). That π + (n)}, we observe that, for 0 ≤ n ≤ N, a ≤ r ≤ b, π + (n, r) is proportional to π is, we have ˜ (n, r, 0), a ≤ r ≤ b, 0 ≤ n ≤ N. π + (n, r) = d π Defining ψ + (n) =

b

π + (n, r), 0 ≤ n ≤ N,

(23)

(24)

r=a

the following lemma gives an expression for the constant d. Lemma 1: The proportionality (or normalizing) constant, d, appearing in (23) is given by ⎤−1

⎡ b n N

a−1 a−1





˜ n−j (−C)−1 e⎦ , π ˜ (n, r, 0)e = g−1 1 − p(n, 0)e = ⎣g + ψ + (j)D d−1 = n=0 r=a

n=0

n=0 j=0

(25) where g = ha

a

ψ + (n)e +

n=0

b

hn ψ + (n)e + hb

n=a+1

N

ψ + (n)e,

(26)

n=b+1

˜ is as given in (2). and D Proof: On dividing (22) by

b N



π ˜ (n, r, 0)e and after some simple algebraic manipulations

n=0 r=a

with the help of (23) and (24), we get b N





π ˜ (n, r, 0)e = g−1 1 −

n=0 r=a

a−1

n=0

11

p(n, 0)e .

(27)

By a similar operation on (4) and using (27), we get

a−1

−1 p(n, 0)e ψ + (0)(−C)−1 . 1− p(0, 0) = g

(28)

n=0

By a similar operation on (5) and using (28), we get ⎤

⎡ n a−1



˜ n−j ⎦ (−C)−1 , 0 ≤ n ≤ a − 1. p(n, 0)e ⎣ ψ + (j)D p(n, 0) = g −1 1 − n=0

(29)

j=0

The stated result follows immediately by post-multiplying (29) by e and using (27). Now we are ready to state the main result of this section. Theorem 1: The steady-state probability vectors {p(n, 0), π(n, r)} and {π + (n, r)} are related by p(n, 0) = E ∗−1

n 

 ˜ n−j (−C)−1 , 0 ≤ n ≤ a − 1, ψ + (j)D

(30)

j=0

    π(0, a) = E ∗−1 π + (0, a) − ψ + (a) − p(a − 1, 0)D C −1 ,    π(0, r) = E ∗−1 π + (0, r) − ψ + (r) C −1 , a + 1 ≤ r ≤ b,   π(n, r) = E ∗−1 π + (n, r) − π(n − 1, r)D C −1 , a ≤ r ≤ b, 1 ≤ n ≤ N − 1,     π(n, b) = E ∗−1 π + (n, b) − ψ + (n + b) − π(n − 1, b)D C −1 , 1 ≤ n ≤ N − b,   π(n, b) = E ∗−1 π + (n, b) − π(n − 1, b)D C −1 , N − b + 1 ≤ n ≤ N − 1, where E∗  g +

n a−1



˜ n−j (−C)−1 e, ψ + (j)D

(31) (32) (33) (34) (35)

(36)

n=0 j=0

and g is as given in (26). Proof From equation (29) and using results of Lemma 1 we have the desired result of (30). Now setting s = 0 in (13) - (17), we get π ˜ (0, a, 0) = π(0, a)C +

b

π ˜ (a, j, 0) + p(a − 1, 0)D,

(37)

π ˜ (r, j, 0), a + 1 ≤ r ≤ b,

(38)

j=a

π ˜ (0, r, 0) = π(0, r)C +

b

j=a

π ˜ (n, r, 0) = π(n, r)C + π(n − 1, r)D, a ≤ r ≤ b − 1, 1 ≤ n ≤ N − 1,

12

(39)

π ˜ (n, b, 0) = π(n, b)C +

b

π ˜ (n + b, j, 0) + π(n − 1, b)D, 1 ≤ n ≤ N − b,

(40)

j=a

π ˜ (n, b, 0) = π(n, b)C + π(n − 1, b)D, N − b + 1 ≤ n ≤ N − 1, Dividing (37) - (41) by

N

b

(41)

˜ π(n, r, 0)e and using (23), (24) and Lemma 1, after little algebraic

n=0 r=a

operations we obtain the desired result (31) - (35). Remarks: (1) Note that Theorem 1 does not have any expression for the probability vectors π(N, r), a ≤ r ≤ b. However, one can obtain the probability vector, π f ull , of the buffer being full at an arbitrary time using the following easily verifiable equation. π f ull =

b

π(N, r) = δ −

r=a

a−1

p(n, 0) −

N −1

b

π(n, r).

(42)

n=0 r=a

n=0

(2) The evaluation of π(N, r)e, a ≤ r ≤ b, is discussed in section 2.2.1. 2.2.1. Evaluation of π(N, r)e, a ≤ r ≤ b As mentioned earlier that the probability vectors, π(N, r), a ≤ r ≤ b, cannot be obtained using the relations given in Theorem 1. However, some probabilities useful in defining key system performance measures can be obtained. Here we will show how one can evaluate 

π(N, r)e, a ≤ r ≤ b. In the following π ∗ (n, r, 0) will denote the derivative of π ∗ (n, r, s) with respect to s evaluated at s = 0. Differentiating the functions in (13-18) with respect to s, setting s = 0 and post-multiplying the resulting ones by e and using (C + D)e = 0 we get ∗

π (0, a, 0)De = π(0, a)e − ha

b

π ˜ (a, j, 0)e − ha p(a − 1, 0)De,

(43)

j=a ∗

π (0, r, 0)De = π(0, r)e − hr

b

π ˜ (r, j, 0)e,

a + 1 ≤ r ≤ b,

(44)

j=a 



π ∗ (n, r, 0)De = π(n, r)e + π ∗ (n − 1, r, 0)De, a ≤ r ≤ b − 1, 1 ≤ n ≤ N − 1, ∗

π (n, b, 0)De = π(n, b)e − hb

b



π ˜ (n + b, j, 0)e + π ∗ (n − 1, b, 0)De, 1 ≤ n ≤ N − b,

(45) (46)

j=a 



π ∗ (n, b, 0)De = π(n, b)e + π ∗ (n − 1, b, 0)De, N − b + 1 ≤ n ≤ N − 1, 

π ∗ (N − 1, r, 0)De + π(N, r)e = 0, a ≤ r ≤ b,

13

(47) (48)



where hr = −Hr∗ (0), a ≤ r ≤ b, is the mean service time for a batch of size r. Using Theorem 1 and the fact that a−1  

π ˜ (n, r, 0)e = g−1 1 − p(n, 0)e π + (n, r)e, a ≤ r ≤ b, 0 ≤ n ≤ N,

(49)

n=0

the equations in (43-48) yield a recursive procedure to solve π(N, r)e, a ≤ r ≤ b, in terms of known quantities. For easy implementation the entire procedure is described below in steps. ˜ Step 1. For a ≤ r ≤ b and 0 ≤ n ≤ N , calculate π(n, r, 0)e using (49). 

Step 2. For a ≤ r ≤ b, calculate π ∗ (0, r, 0)De, using (43) and (44). Step 3. For n from 1 to N − 1, and for r from a to b − 1 

calculate π ∗ (n, r, 0)De using (45). 

Step 4. For n from 1 to N − b calculate π ∗ (n, b, 0)De using (46). 

Step 5. For n from N − b + 1 to N − 1 calculate π ∗ (n, b, 0)De using (47). Step 6. Finally, calculate π(N, r)e, a ≤ r ≤ b, using (48). 2.3. Joint distribution at an arrival epoch In this section we will give an expression for the steady-state probability vector at an arrival epoch, i.e., joint distribution of the number of customers in the queue, the number in service, and the phase of the arrival process at an arrival epoch. Suppose that the ith component of p− (n, 0), 0 ≤ n ≤ a − 1, gives the steady-state probability that an arrival finds n, 0 ≤ n ≤ a − 1, customers in the queue with the server being idle and the arrival process is in state i, and that the ith component of π − (n, r), a ≤ r ≤ b, 0 ≤ n ≤ N, gives the steady-state probability that an arrival finds n, 0 ≤ n ≤ N, customers in the queue with the server busy in serving a batch of r, a ≤ r ≤ b, customers and the arrival process is in state i. Then the vectors p− (n, 0) and π − (n, r) are given by (see, Gupta and Laxmi [28]) p− (n, 0) =

1 λ∗ p(n, 0)D,

0 ≤ n ≤ a − 1, (50)

π − (n, r)

=

1 λ∗ π(n, r)D,

a ≤ r ≤ b, 0 ≤ n ≤ N − 1.

Remarks: (1) Note that it is possible for an arrival to be lost when the buffer is full. Thus, the effective arrival rate (λe ) is given by λe = [λ∗ −

b

r=a

14

π(N, r)De],

(51)

which is equal to the effective service rate in steady state. (2) Since π(N, r), a ≤ r ≤ b, is not known we cannot compute π − (N, r), a ≤ r ≤ b. However, b  π − (N, r) which is we can compute quantities such as an arrival is lost, as we need only obtained as

1 f ull D λ∗ π

r=a

using (42).

3. Performance measures In this section we will list a number of system performance measures of interest in steady state along with their expressions. These are in addition to the ones mentioned in earlier sections. 1. The probability mass function of the number of customers in the system at an arbitrary epoch is given by

(n)

psystem =

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩

p(n, 0)e, min(b,n)  r=max(n−N,a)

0 ≤ n ≤ a − 1,

π(n − r, r)e, a ≤ n ≤ N + b.

2. The probability mass function of the number of customers in the queue at an arbitrary epoch is given by

p(n) queue =

⎧ b  ⎪ ⎪ ⎪ p(n, 0)e + π(n, r)e, ⎪ ⎪ ⎨ r=a ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

b  r=a

0 ≤ n ≤ a − 1,

π(n, r)e, a ≤ n ≤ N.

3. The (conditional) probability mass function of the number of customers in service given that the server is busy is (r) pservice

=c

N

π(n, r)e, a ≤ r ≤ b,

n=0

where c is the normalizing condition given by c−1 = [1 −

a−1  n=0

p(n, 0)e] = Pbusy

4. The probability, Pbusy , that the server is busy at an arbitrary epoch is given by Pbusy =

N

b

π(n, r)e.

n=0 r=a

5. The probability, Plost , that an arriving customer will be lost is given by Plost =

b 1

π(N, r)De. λ∗ r=a

15

6. The mean, L, number of customers in the system can be obtained using the probability (n)

function psystem . 7. The mean, Lq , number of customers in the queue can be obtained using the probability (n)

function pqueue . 8. The mean, Lservice , number of customers serviced in a batch can be obtained using the (r)

probability function pservice . 9. Using Little’s law we can obtain the mean, W , waiting time in the system as well as the mean, Wq , waiting time in the queue as W =

L λe

and Wq =

Lq λe ,

where λe is the effective

arrival rate as given in (51). 4. Phase type services In the case when the service time is modeled as a phase type distribution (PH-distribution), the computations become much simpler and use only matrix arithmetic. Recall that a PHdistribution is obtained as the time until absorption in a finite state Markov chain with one absorption state. It is characterized by an initial probability vector and a square matrix governing the transitions to various transient states. PH-distributions are defined for both discrete and continuous time. For details on PH-distributions and their properties, we refer the reader to Neuts [34, 43]. The following theorem shows the computational aspects in the case of PHservices. Theorem 2: Suppose that Hr (.), a ≤ r ≤ b, follows a PH-distribution with an irreducible representation (βr , Sr ) of dimension ν. [Note that we assume that all PH-distributions are of the same dimension ν but can be modified to include other cases and the analysis is very similar]. (r) (b) (r) ˜n,N +1 appearing in (2), are given by , Bn,k , and B Then the matrices A , A˜ , A˜ k

k

N +1

(r)

Ak = (Im ⊗ βr )Mk,r (Im ⊗ Sr0 ), 0 ≤ k ≤ N, a ≤ r ≤ b,   −1  (b) (Im ⊗ Sb0 ), b ≤ k ≤ N, A˜k = (Im ⊗ βb ) − Mk−1,b (D ⊗ Iν ) (C ⊕ Sb ) + (D ⊗ Iν )   −1  (r) (Im ⊗ Sr0 ), a ≤ r ≤ b, A˜N +1 = (Im ⊗ βr ) − MN,r (D ⊗ Iν ) (C ⊕ Sr ) + (D ⊗ Iν ) (a)

˜ a−n A , 0 ≤ n ≤ a − 1, 0 ≤ k ≤ N, Bn,k = D k ˜ a−n A˜(a) , 0 ≤ n ≤ a − 1 , ˜n,N +1 = D B N +1

16

(52)

where

Sr0 = −Sr e, Mk,r = Mk−1,r (D ⊗ Iν )M0,r , 1 ≤ k ≤ N, a ≤ r ≤ b,

(53)

M0,r = −(C ⊕ Sr )−1 , a ≤ r ≤ b, ˜ is as given in (2). and D Proof: First note that the matrices {P˜ (n, t), n ≥ 0, t ≥ 0} associated with the counting process for M AP satisfy (see, e.g., Neuts and Li [47]) P˜  (n, t) = P˜ (n, t)C + P˜ (n − 1, t)D, n ≥ 0,

(54)

with P˜ (−1, t) = 0 and P˜ (0, 0) = I. Using the properties of Kronecker product and noticing that P˜ (k, t) → 0 as t → ∞, we see that (r) Ak

∞ =

P˜ (k, t)dHr (t) = (Im ⊗ βr )Mk,r (Im ⊗ Sr0 ), 0 ≤ k ≤ N, a ≤ r ≤ b,

(55)

0

where Mk,r =

∞ 0

=

P˜ (k, t) ⊗ eSr t dt   ⎧  −1 − ∞ P ˜  (0, t) ⊗ eSr t S −1 dt, k = 0, a ≤ r ≤ b, ⎪ ⊗ S − I ⎪ m r r 0 ⎨ ⎪ ⎪ ⎩

−

∞ 0

(56)

P˜  (k, t) ⊗ eSr t Sr−1 dt, 1 ≤ k ≤ N, a ≤ r ≤ b.

Using (54) in (56) and after some routine manipulations involving Kronecker products and sums, we get

M0,r = −(C ⊕ Sr )−1 , a ≤ r ≤ b, (57) Mk,r = Mk−1,r (D ⊗ Iν )M0,r , 1 ≤ k ≤ N, a ≤ r ≤ b.

On noticing that (b) A˜k =

where

∞ 

∞

j=k

(b)

˜ b (Im ⊗ S 0 ), b ≤ k ≤ N, Aj = (Im ⊗ βb )M b

(58)

∞

˜b = M

P˜ (j, t) ⊗ eSb t dt = −Mk−1,b (D ⊗ Iν )[(C ⊕ Sb ) + (D ⊗ Iν )]−1 ,

j=k 0

17

(59)

(b)

stated expression for A˜k holds good. By a similar calculation we get the stated expression (r) for A˜N +1 , a ≤ r ≤ b. From these the expressions for Bn,k , 0 ≤ n ≤ a − 1, 0 ≤ k ≤ N, and ˜n,N +1 , 0 ≤ n ≤ a − 1, follow immediately. The detailed proof of Theorem 2 is given in B Appendix B. Note:- In this section, we briefly discussed the necessary steps required for the computation of (r)

the elements, the matrices Ak , Bn,k , etc., of the TPM (P ). The evaluation of these matrices, in general, for arbitrary service time distribution requires numerical integration. However, for phase type service time distribution, it is not necessary to compute integrals. After specifying the TPM completely, for PH service time distribution, one has to solve (3) to obtain the departureepoch distributions. 5. Numerical results The purpose of this section is to bring out the qualitative aspects of the queueing system under consideration through illustrative numerical examples. For the arrival process, we consider the following five sets of values for C and D. 1. Erlang (ERA): Erlang of order 2.     −2 2 0 0 C= ,D = 0 −2 2 0 2. Exponential (EXA): C=



−1



,D =



1



3. Hyperexponential (HEA):     −1.90 0 1.710 0.190 C= ,D = 0 −0.19 0.171 0.019 4. M AP with negative correlation (M N A): ⎞ ⎛ ⎞ ⎛ 0 0 0 −1.00222 1.00222 0 ⎟ ⎜ ⎟ ⎜ ⎟ , D = ⎜ 0.01002 0 0.99220 ⎟ C=⎜ 0 −1.00222 0 ⎠ ⎝ ⎠ ⎝ 223.49250 0 2.25750 0 0 −225.75 5. M AP with positive correlation (M P A): ⎛ −1.00222 1.00222 0 ⎜ ⎜ C=⎝ 0 −1.00222 0 0

0

⎞

⎛

0

0

0

⎞

⎟ ⎜ ⎟ ⎟ , D = ⎜ 0.99220 0 0.01002 ⎟ ⎠ ⎝ ⎠ −225.75 2.25750 0 223.49250

All these five M AP processes will be normalized so as to have a specific arrival rate. Note that these are qualitatively different in that they have different variance and correlation structure. 18

The first three arrival processes, namely ERA, EXA, and HEA, correspond to renewal processes and so the correlation is 0. The arrival processes labelled M N A has correlated arrivals with correlation between two successive inter-arrival times given by -0.4889, and the arrivals corresponding to the process labelled M P A has positive correlation with a value 0.4889. The ratios of the standard deviations of the inter-arrival times of these five arrival processes with respect to ERA are, respectively, 1, 1.4142, 3.1745, 1.9934, and 1.9934. For the service times (βr , Sr ), we consider the following three P H−distributions. A. Erlang (ERS):

 βr = (1, 0), Sr =

B. Exponential (EXS): β r = 1, Sr =



−2μr

2μr

0

−2μr



−μr



.

.

C. Hyperexponential (HES):  β r = (0.9, 0.1), Sr =



−1.9μr

0

0

−0.19μr

.

These PH-representations are such that they all have the same mean of of r, a ≤ r ≤ b, customers. However, they are qualitatively different.

1 μr

when serving a batch

We consider the following three scenarios for the service rates, μr , a ≤ r ≤ b. I. Constant (CON ): Here we take μr = μ, a ≤ r ≤ b. Thus, the service rates are independent of the batch size. II. Inversely proportional (IP ): Here we take μr =

μ r,a

≤ r ≤ b, so that the process-

ing time will increase as the batch size increases. This may be the case in some applications. III. Linearly decreasing (LD): Here we take μr = μ1 −

μ1 −μ2 b−a (r

− a), a ≤ r ≤ b, where

μ1 ≥ μ2 > 0, so that the service rate comprises of a fixed rate and a rate that decreases as the batch size increases. This again may be appropriate for some applications. In order to compare various scenarios we need to normalize appropriate parameters so as to have the same weighted average service time which is reciprocal of the weighted average service  rate and it is calculated using the formula

b

rμr . r=a r

r=a  b

19

For the rest of this section, we fix λ∗ = 1 and take the service rates for the three scenarios, I, II, and III, as given in Table 1. Table 1: Service rates for the three scenarios Batch size(r)

μr

I

II

III

5

μ5

0.5

0.728571

1.273810

6

μ6

0.5

0.607143

1.019048

7

μ7

0.5

0.520408

0.764286

8

μ8

0.5

0.455357

0.509524

9

μ9

0.5

0.404762

0.254762

Note that all three scenarios yield an average service time of a batch to be 2 and the buffer size is considered to be 15. We look at the impact of the arrival process, the service time distributions, and the three scenarios for the service rates on some selected system performance measures, which are displayed in Table 2. Table 2: Selected measures for various combinations of arrivals and services ERS M easures M AP Plost

I

II

EXS III

I

II

HES III

I

II

III

ERA 0.0000 0.0000 0.0000 0.0005 0.0000 0.0000 0.0474 0.0215 0.0064 EXA 0.0001 0.0000 0.0000 0.0009 0.0001 0.0000 0.0497 0.0234 0.0076 HEA 0.0012 0.0003 0.0001 0.0060 0.0020 0.0007 0.0649 0.0364 0.0156 MNA 0.0001 0.0000 0.0000 0.0010 0.0002 0.0000 0.0499 0.0235 0.0077 MPA 0.4230 0.4187 0.4139 0.4202 0.4155 0.4112 0.4226 0.4143 0.4094

Lq

ERA 2.0815 2.0152 2.0007 2.2213 2.0614 2.0052 3.1816 2.6742 2.2587 EXA 2.1207 2.0313 2.0029 2.2694 2.0896 2.0121 3.1972 2.6936 2.2804 HEA 2.2515 2.0932 2.0144 2.4477 2.2051 2.0501 3.2572 2.7843 2.3854 MNA 2.1530 2.0458 2.0059 2.3043 2.1091 2.0180 3.2165 2.7075 2.2894 MPA 2.1563 2.1437 2.1988 2.1875 2.1675 2.2306 2.5292 2.3553 2.3335

Lservice

ERA 5.0741 5.0197 5.0008 5.1680 5.0786 5.0098 5.2698 5.3128 5.3532 EXA 5.1305 5.0560 5.0077 5.2229 5.1330 5.0335 5.2889 5.3310 5.3899 HEA 5.3903 5.2656 5.0797 5.4732 5.3990 5.2073 5.3891 5.4370 5.5804 MNA 5.1421 5.0642 5.0103 5.2341 5.1435 5.0398 5.2924 5.3343 5.3946 MPA 5.2681 5.4547 6.0987 5.2893 5.4761 6.1242 5.3876 5.5737 6.1687

An examination of the table reveals a number of observations with respect to the measures, Plost , Lq , and Lservice . First looking at Plost , we notice the following. • While this measure is significantly small for all but M P A arrivals in the case of Erlang (ERS) and exponential (EXS) services, we notice this to be significant for the hyperexponential (HES) services. In the case of M P A arrivals, the measure is significantly large and also appears to be insensitive to the type of service distribution used for the services. 20

• As we change the scenario we notice that this measure appears to decrease from Scenario I to II to III. This phenomenon may look counter intuitive initially since having a constant service rate (as opposed to the rates decreasing as the batch size increases) one will expect to see a smaller probability for lost customers. However, to preserve the same weighted average service time, the rates are initially higher for Scenarios II and III as compared to Scenario I. This explains the phenomenon. Also, the rate of decrease (as a percentage) is not significant in the case of M P A arrivals, while we see this to be very significant for HEA arrivals. In the case Lq we observe that this measure appears to show a similar behavior for all arrivals and for the three scenarios in the case of ERS and EXS services. However, for the HES services, we see that all but M P A arrivals this measure has a higher value (as compared to the other two services) and for M P A arrivals there is only a slight increase. With regard to Lservice we notice the following. • For all three scenarios and for all services, we notice that this measure appears to be closer to the minimum batch size, namely, 5, in the case of ERA, EXA, HEA, and M N A arrivals. In the case of M P A arrivals this measure is greater than 6 for Scenario III; however for the other two scenarios it is less than 6 and is comparable to HEA arrivals. • For all arrivals, we see an interesting pattern. By first looking at all but M P A arrivals, we see that this measure appears to decrease as we go from Scenario I to II to III when the services are either ERS or EXS type. However, for the HES services, the trend is reversed. That is, for this case, this measure appears to increase as we move from Scenario I to II to III. If we now look at M P A arrivals, we see that this measure appears to increase when moving from Scenario I to II to III. This appears to be the case for all services. The observations recorded above show not only the significant role of the variability in the arrivals and services but also on the correlation, especially, the positive one. 6. Concluding remarks In this paper, we studied a finite capacity bulk service queue with Markovian arrival process. The service times, which depend on the size of the batches, are generally distributed. When the services are of phase type we developed efficient procedures for computing the steady-state probabilities at various epochs. An illustrative example to show the impact of the correlated arrivals on selected system performance measures is presented. The procedure developed in this paper can be used to analyze more complex queuing models involving BM AP arrivals and/or batch Markovian service process (BM SP ).

21

Acknowledgement The authors thank the anonymous referee for their valuable comments and suggestions that helped to improve the presentation of the paper in the current form. The second author acknowledge the Department of Science and Technology (DST), Govt. of India, for the financial support under the project grant SR/S4/MS:789/12. [1] L. Abolnikov, A. Dukhovny, Optimization in HIV screening problems, Journal of Applied Mathematics and Stochastic Analysis 16 (4) (2003) 361–374. [2] V. Goswami, J. Mohanty, S. K. Samanta, Discrete-time bulk-service queues with accessible and non-accessible batches, Applied Mathematics and Computation 182 (1) (2006) 898 – 906. [3] S. K. Bar-Lev, M. Parlar, D. Perry, W. Stadje, F. A. V. der Duyn Schouten, Applications of bulk queues to group testing models with incomplete identification, European Journal of Operational Research 183 (1) (2007) 226–237. [4] D. Claeys, J. Walraevens, K. Laevens, H. Bruneel, A queueing model for general group screening policies and dynamic item arrivals, European Journal of Operational Research 207 (2) (2010) 827–835. [5] D. Claeys, J. Walraevens, K. Laevens, H. Bruneel, Analysis of threshold-based batch-service queueing systems with batch arrivals and general service times, Performance Evaluation 68 (6) (2011) 528–549. [6] D. Claeys, B. Steyaert, J. Walraevens, K. Laevens, H. Bruneel, Tail distribution of the delay in a batch-service queueing model, Computers & Operations Research 39 (11) (2012) 2733–2741. [7] I. J. B. F. Adan, J. A. C. Resing, Multi-server batch-service systems, Statistics Neerlandica 54 (2) (2000) 202–220. [8] R. Arumuganathan, S. Jeyakumar, Steady state analysis of a bulk queue with multiple vacations, setup times with N-policy and closedown times, Applied Mathematical Modelling 29 (10) (2005) 972–893. [9] Y. Baba, A bulk service GI/M/1 queue with service rates depending on service batch size, Journal of Operations Research 39 (1) (1996) 25–34. [10] A. D. Banik, M. L. Chaudhry, U. C. Gupta, On the finite buffer queue with renewal input and batch Markovian service process: GI/BM SP/1/N , Methodology and Computing in Applied Probability 10 (8) (2008) 559–575.

22

[11] A. D. Banik, U. C. Gupta, M. L. Chaudhry, Finite-buffer bulk service queue under Markovian service process: GI/M SP (a,b) /1/N , Journal of Stochastic Analysis and Application 27 (3) (2009) 500–522. [12] M. L. Chaudhry, J. G. C. Templeton, A First Course in Bulk Queues, John Wiley & Sons., New York, 1983. [13] M. L. Chaudhry, U. C. Gupta, Modelling and analysis of M/G(a,b) /1/N queue - A simple alternative approach, Queueing Systems 31 (3) (1999) 95–100. [14] H. Gold, P. Tran-Gia, Performance analysis of a batch service queue arising out of manufacturing system modeling, Queueing Systems 14 (3) (1993) 413–426. [15] G. H´ebuterne, C. Rosenberg, Arrival and departure state distributions in the general bulkservice queue, Naval Research Logistics 46 (1) (1999) 107–118. [16] M. F. Neuts, A general class of bulk queues with Poisson input, Annals of Mathematical Statistics 38 (3) (1967) 759–770. [17] W. B. Powell, P. Humblet, The bulk service queue with a general control strategy: Theoretical analysis and a new computational procedure, Operations Research 34 (2) (1986) 267–275. [18] K. Sikdar, U. C. Gupta, Analytic and numerical aspects of batch service queues with single vacation, Computers & Operations Research 32 (4) (2005) 943–966. [19] G. L. Curry, R. M. Feldman, An M/M/1 queue with a general bulk service rule, Naval Research Logistics Quarterly 32 (3) (1985) 595–603. [20] M. F. Neuts, Transform-free equations for the stationary waiting time distributions in the queue with Poisson arrivals and bulk services, Annals of Operations Research 8 (3) (1987) 3–26. (a,b)

[21] M. L. Chaudhry, J. Gai, A simple and extended computational analysis of M/Gj (a,b) M/Gj /1/B

/1 and

+ b queues using roots, INFOR 50 (2) (2012) 72–79.

[22] A. Banerjee, U. C. Gupta, Reducing congestion in bulk-service finite-buffer queueing system using batch-size-dependent service, Performance Evaluation 69 (1) (2012) 53–70. [23] R. Germs, N. van Foreest, Analysis of finite-buffer state-dependent bulk queues, OR Spectrum 35 (3) (2013) 563–583. [24] D. Claeys, B. Steyaert, J. Walraevens, K. Laevens, H. Bruneel, Tail probabilites of the delay in a batch-service queueing model with batch-size dependent service times and a timer mechanism, Computers & Operations Research 40 (2013) 1497–1505. 23

[25] F. Avram, A. G´omez-Corral, On bulk-service M AP/P H L,N /1/N G-Queues with repeated attempts, Annals of Operations Research 141 (1) (2006) 109–137. [26] A. D. Banik, Queueing analysis and optimal control of BM AP/G(a,b) /1/N and BM AP/M SP (a,b) /1/N systems, Computers & Industrial Engineering 57 (3) (2009) 748– 761. [27] S. Chakravarthy, Analysis of a finite M AP/G/1 queue with group services, Queueing Systems 13 (4) (1993) 385–407. [28] U. C. Gupta, P. V. Laxmi, Analysis of the M AP /Ga,b /1/N queue, Queueing Systems 38 (2) (2001) 109–124. [29] M. L. Chaudhry, U. C. Gupta, Analysis of a finite-buffer bulk-service queue with discretemarkovian arrival process: D − M AP /G(a,b) /1/N , Neval Research Logistics 50 (4) (2003) 345–363. [30] K. Sikdar, U. C. Gupta, The queue length distributions in the finite buffer bulk-service M AP /G/1 queue with multiple vacation, Sociedad de Estadistica e Investigacidn Operativa Top 13 (1) (2005) 75–103. [31] K. Sikdar, Analysis of the M AP /G(a,b) /1/N queue with multiple vacations, Applied Mathematical Modelling 32 (2008) 1308–1317. [32] D. Claeys, B. Steyaert, J. Walraevens, K. Laevens, H. Bruneel, Analysis of a versatile batchservice queueing model with correlation in the arrival process, Performance Evaluation 70 (4) (2013) 300–316. [33] D. M. Lucantoni, New results on the single server queue with a batch markovan arrival process, Stochastic Models 7 (1) (1991) 1–46. [34] M. F. Neuts, Structured Stochastic Matrices of M/G/1 Type and Their Applications, Marcel Dekker Inc., New York, 1989. [35] S. R. Chakravarthy, The batch Markovian arrival process: A review and future work, Advances in Probability Theory and Stochastic Processes. Eds., A. Krishnamoorthy et al., Notable Publications Inc., NJ, 2001. [36] A. D. Banik, U. C. Gupta, S. S. Pathak, BMAP/G/1/N queue with vacations and limited service discipline, Applied Mathematics and Computation 180 (2) (2006) 707 – 721. [37] U. C. Gupta, K. Sikdar, Computing queue length distributions in MAP/G/1/N queue under single and multiple vacation, Applied Mathematics and Computation 174 (2) (2006) 1498–1525. 24

[38] S. R. Chakravarthy, Markovian Arrival Processes, Wiley Encyclopedia of Operations Research and Management Science, Published Online: 15 JUN 2010, 2010. [39] A. G.-C. J. Artalejo, Q. He, Markovian arrivals in stochastic modelling: a survey and some new results, SORT 34 (2010) 101–144. [40] S. Chakravarthy, A. S. Alfa, A multi-server queue with markovian arrivals and group services with thresholds, Naval Research Logistics 40 (1993) 811–827. [41] S. Chakravarthy, A. S. Alfa, A finite capacity queue with Markovian arrivals and two servers with group services, Journal of Applied Mathematics and Stochastic Analysis 7 (1994) 161– 178. [42] A. S. Alfa, K. L. Dolhun, S. R. Chakravarthy, A discrete single server queue with Markovian arrivals and phase type group services, Journal of Applied Mathematics and Stochastic Analysis 8 (1995) 151–176. [43] M. F. Neuts, Matrix-geometric solutions in stochastic models : An algorithmic approach, Johns Hopkins University Press, Baltimore, 1981. [44] S. R. Chakravarthy, L. Bin, A finite capacity queue with nonrenewal input and exponential dynamic group services, ORSA Journal on Computing 9 (1997) 276–287. [45] M. Maecus, H. Minc, A survey of matrix theory and matrix inequalities, Allyn and Bacon, Boston, MA, 1964. [46] D. M. Lucantoni, K. S. Meier-Hellstern, M. F. Neuts, A single-server queue with server vacations and a class of non-renewal arrival processes, Advances in Applied Probability 22 (3) (1990) 676–705. [47] M. F. Neuts, J.-M. Li, An algorithm for the P (n, t) matrices of a continuous BMAP, In Matrix-analytic Methods in Stochastic Models, Eds., S.R. Chakravarthy and A.S. Alfa, Marcel Dekker, 1996.

25

Appendix A. Derivation of equation (4) to (11) Relating the states of the system at time t and t + dt and using the supplementary variable technique, we have for 1 ≤ i ≤ m and u ≥ 0 m b



d pˆi (0, 0; t) = pˆk (0, 0; t)cki + π ˆi (0, j, 0; t), dt j=a

k=1

d pˆi (n, 0; t) = dt

m

(A.1)

pˆk (n, 0; t)cki +

k=1

m

pˆk (n − 1, 0; t)dki +

b

π ˆi (n, j, 0; t),

j=a

k=1

1 ≤ n ≤ a − 1,

m m



∂  ˜ a (u) − π ˆi (0, a, u; t) = π ˆk (0, a, u; t)cki + pˆk (a − 1, 0, t)dki h ∂t ∂u

∂

k=1

k=1

b

+

(A.2)

˜ a (u), π ˆi (a, j, 0; t)h

(A.3)

j=a

∂ ∂t ∂ ∂t ∂ ∂t

− −

−

m b



∂  ˜ r (u), π ˆi (0, r, u; t) = π ˆk (0, r, u; t)cki + π ˆi (r, j, 0; t)h ∂u

∂  π ˆi (n, r, u; t) = ∂u ∂  π ˆi (n, b, u; t) = ∂u

k=1 m

π ˆk (n, r, u; t)cki +

j=a m

k=1

k=1

m

m

π ˆk (n, b, u; t)cki +

k=1

a + 1 ≤ r ≤ b,

π ˆk (n − 1, r, u; t)dki ,

(A.5)

a ≤ r ≤ b − 1, 1 ≤ n ≤ N − 1, b

˜ b (u), π ˆk (n − 1, b, u; t)dki + π ˆi (n + b, j, 0; t)h j=a

k=1

1 ≤ n ≤ N − b,

m m



∂  − π ˆi (n, b, u; t) = π ˆk (n, b, u; t)cki + π ˆk (n − 1, b, u; t)dki ∂t ∂u

∂

k=1

(A.4)

(A.6) (A.7)

k=1

N − b + 1 ≤ n ≤ N − 1, m m m 



∂ − π ˆi (N, r, u; t) = π ˆk (N, r, u; t)cki + π ˆk (N, r, u; t)dki + π ˆk (N − 1, r, u; t)dki , ∂t ∂u

∂

k=1

k=1

k=1

a ≤ r ≤ b. Hence in steady-state, letting t → ∞, and using the fact that pi (n, 0) = lim pˆi (n, 0; t), t→∞

ˆi (n, r, u; t), π ˜i (n, r, u) = lim π t→∞

1 ≤ i ≤ m, 0 ≤ n ≤ a − 1, 1 ≤ i ≤ m, 0 ≤ n ≤ N, a ≤ r ≤ b.

26

(A.8)

the equations (A.1) - (A.8) reduces to 0 = 0 =

m

pk (0, 0)cki +

b

k=1

j=a

m

m

pk (n, 0)cki +

k=1

π ˜i (0, j, 0),

(A.9)

pk (n − 1, 0)dki +

b

π ˜i (n, j, 0), 1 ≤ n ≤ a − 1,

j=a

k=1

(A.10) − − −

d π ˜i (0, a, u) = du d π ˜i (0, r, u) = du d π ˜i (n, r, u) = du

m

π ˜k (0, a, u)cki +

m

k=1

k=1

m

b

π ˜k (0, r, u)cki +

k=1 m

π ˜k (n, r, u)cki +

k=1

˜ a (u) + pk (a − 1, 0)dki h

b

˜ a (u),(A.11) π ˜i (a, j, 0)h

j=a

˜ r (u), π ˜i (r, j, 0)h

j=a m

a + 1 ≤ r ≤ b,

π ˜k (n − 1, r, u)dki ,

k=1

a ≤ r ≤ b − 1, 1 ≤ n ≤ N − 1, −

− −

(A.12)

(A.13)

m m b



d ˜ b (u), π ˜i (n, b, u) = π ˜k (n, b, u)cki + π ˜k (n − 1, b, u)dki + π ˜i (n + b, j, 0)h du

d π ˜i (n, b, u) = du

d π ˜i (N, r, u) = du

k=1

k=1

m

m

π ˜k (n, b, u)cki +

k=1

m

j=a

1 ≤ n ≤ N − b,

π ˜k (n − 1, b, u)dki , N − b + 1 ≤ n ≤ N − 1,(A.15)

k=1 m

π ˜k (N, r, u)cki +

k=1

(A.14)

π ˜k (N, r, u)dki +

k=1

m

π ˜k (N − 1, r, u)dki ,

k=1

a ≤ r ≤ b.

(A.16)

Then from (A.9) - (A.16), after using matrix and vector notations, we have the set of (matrix) differential equations (4) - (11). Appendix B. Detailed proof of Theorem 2 For 0 ≤ k ≤ N, a ≤ r ≤ b, we have ∞ ∞ ∞ (r) (Sr t) 0 ˜ ˜ = Sr dt = Im P˜ (k, t)Im ⊗β r e(Sr t) S0r dt, P (k, t)dHr (t) = P (k, t)⊗β r e Ak 0

= ∴

(r) Ak

=

 

Im ⊗β r

  ∞

0



P˜ (k, t)⊗e(Sr t) dt

0

  Im ⊗S0r ,

0

   Im ⊗β r Mk,r Im ⊗S0r ,

where

∞ Mk,r =

(B.1)

P˜ (k, t)⊗e(Sr t) dt,

0

27

0 ≤ k ≤ N, a ≤ r ≤ b

(B.2)

Now in order to calculate Mk,r , we integrate RHS of (B.2) by parts, and obtain Mk,r

    ∞ −1 ˜ = − P (k, 0)⊗Sr − P˜  (k, t)⊗e(Sr t) Sr−1 dt, as t → ∞, P˜ (k, t) → 0

(B.3)

0

Now setting k = 0 in (B.3), we obtain for a ≤ r ≤ b M0,r

  ∞ −1 = − P˜ (0, 0)⊗Sr − P˜  (0, t)⊗e(Sr t) Sr−1 dt, 0

 ∞  −1 − P˜ (0, t)C⊗e(Sr t) Sr−1 dt, as P˜  (0, t) = P˜ (0, t)C = − Im ⊗Sr 0

  = − Im ⊗Sr−1 −

∞ 

P˜ (0, t)⊗e(Sr t) dt

0

  ! = − Im ⊗Sr−1 − M0,r C⊗Sr−1 ,

  C⊗Sr−1 ,

  using (B.2) .

Therefore, M0,r + M0,r C⊗Sr−1

!

 −1   , as (A⊗B)−1 = A−1 ⊗B −1 . = − Im ⊗Sr

So the above expression can be now written in the following form   −1 ! −1 . = − Im ⊗Sr M0,r Imν + C⊗Sr 

! Imν + C⊗Sr−1 is invertible. Hence, post   ! −1 we have for a ≤ r ≤ b multiplying both side of the above expression by Imν + C⊗Sr−1

Since both C and Sr are invertible, the term

∴

 −1  !−1 Imν + C⊗Sr−1 M0,r = − Im ⊗Sr  −1 ! = − Imν + C⊗Sr−1 , Im ⊗Sr    −1  ! = − Imν Im ⊗Sr + C⊗Sr−1 Im ⊗Sr   −1 = − Im ⊗Sr + C Im ⊗Sr−1 Sr ,  −1  , = − Im ⊗Sr + C⊗Iν −1  , = − C⊗Iν + Im ⊗Sr   !−1 , as A⊕B = A⊗Im + In ⊗B M0,r = − C ⊕ Sr

28

(B.4)

So (B.4) gives the value of M0,r (a ≤ r ≤ b). Now from (B.3) for 1 ≤ k ≤ N, a ≤ r ≤ b we have Mk,r

  ∞ −1 ˜ = − P (k, 0)⊗Sr − P˜  (k, t)⊗e(Sr t) Sr−1 dt, 0

 ∞    −1 − P˜ (k, t)C + P˜ (k − 1, t)D ⊗e(Sr t) Sr−1 dt, = − 0m ⊗Sr 

0

 ∵ P˜ (k, 0) = 0m , & P˜  (k, t) = P˜ (k, t)C + P˜ (k − 1, t)D, ∀ k ≥ 1 ∞

= −0mν −  ∞ = −

P˜ (k, t)C⊗e(Sr t) Sr−1 dt

0

 (Sr t)

P˜ (k, t)⊗e

dt

∞ − 0



C⊗Sr−1



0

P˜ (k − 1, t)D⊗e(Sr t) Sr−1 dt,

 ∞

−

    = −Mk,r C⊗Sr−1 − Mk−1,r D⊗Sr−1 ,

0



 (Sr t)

P˜ (k − 1, t)⊗e

dt

  D⊗Sr−1 ,

 using (B.2)

The above expression can be written in the following form      = −Mk−1,r D⊗Sr−1 . Mk,r Imν + C⊗Sr−1 −1  we have Now post multiplying both side of the above expression by Imν + C⊗Sr−1

∴

  −1 , Mk,r = −Mk−1,r D⊗Sr−1 Imν + C⊗Sr−1   −1 , = −Mk−1,r Imν + C⊗Sr−1 D −1 ⊗Sr   −1 , = −Mk−1,r D −1 ⊗Sr + C D −1 ⊗Iν      −1 , = −Mk−1,r Im ⊗Sr D −1 ⊗Iν + C⊗Iν D −1 ⊗Iν     −1 , Im ⊗Sr + C⊗Iν D −1 ⊗Iν = −Mk−1,r   −1 , = −Mk−1,r C⊕Sr D −1 ⊗Iν   −1 −1 , = −Mk−1,r C⊕Sr D⊗Iν   −1  , = −Mk−1,r D⊗Iν C⊕Sr    −1  − C⊕Sr , = Mk−1,r D⊗Iν   ! using (B.4) . Mk,r = Mk−1,r D⊗Iν M0,r ,

29

(B.5)

Since M0,r (a ≤ r ≤ b) are known from (B.4), hence the values of Mk,r (1 ≤ k ≤ N, a ≤ r ≤ b) can be calculated recursively from (B.5). Now (b) A˜k =

= ˜b = Where M

∞

j=k





˜b = M

(b)

Aj =

Im ⊗β b

∞ j=k

∞ 



∞  ∞     

 

Mj,b Im ⊗Sb0 Im ⊗β b Mj,b Im ⊗Sb0 = Im ⊗β b j=k

  0 ˜ Mb Im ⊗Sb ,

j=k

b ≤ k ≤ N.

(B.6)

 Mj,b , which we have to evaluate. Hence, from (B.6), we have Mj,b



∞     

= Mj−1,b D⊗Iν M0,b , using (B.5)

j=k

j=k

= −

∞ 

Mj−1,b

  −1 , D⊗Iν C ⊕ Sb

  using (B.4)

j=k

∞   −1 

Mj−1,b D⊗Iν C ⊕ Sb , = − Mk−1,b + j=k+1

  −1  ˜ b D⊗Iν C ⊕ Sb . = − Mk−1,b + M

∴

   −1    −1 ˜ b Imν + D⊗Iν C ⊕ Sb . = −Mk−1,b D⊗Iν C ⊕ Sb M

  −1 −1  , we Now post multiplying both side of the above expression by Imν + D⊗Iν C ⊕ Sb obtain    −1   −1 −1 ˜ b = −Mk−1,b D⊗Iν C ⊕ Sb , Imν + D⊗Iν C ⊕ Sb M     −1  −1 , Imν + D⊗Iν C ⊕ Sb C ⊕ Sb = −Mk−1,b D⊗Iν     −1 . C ⊕ Sb + D⊗Iν = −Mk−1,b D⊗Iν Therefore from (B.6) we obtain (b) A˜k =

      −1   Im ⊗β b − Mk−1,b D⊗Iν C ⊕ Sb + D⊗Iν Im ⊗S0b , b ≤ k ≤ N

 (r) (r) Similarly, one can obtain the results for A˜N +1 = ∞ j=N +1 Aj (a ≤ r ≤ b). (a) (a) ˜n,k and B ˜n,N +1 can be obtained as are known now, B Since A , A˜ k

N +1

˜ a−n A(a) , 0 ≤ n ≤ a − 1, 0 ≤ k ≤ N, ˜n,k = D B k ˜ a−n A˜(a) , 0 ≤ n ≤ a − 1 . ˜n,N +1 = D B N +1

30