# Computational Procedures for a Class of GI/D/k Systems in Discrete Time

Authors: Md. Mostafizur Rahman; Attahiru Sule Alfa
Journal: Journal of Probability and Statistics
 (2009)

---

## Abstract

A class of discrete time GI/D/k
systems is considered for which the interarrival times have finite support and customers are served in first-in first-out (FIFO) order. The system is formulated as a single server queue with new general independent interarrival times and constant service duration by assuming cyclic assignment of customers to the identical servers. Then the queue length is set up as a quasi-birth-death (QBD) type Markov chain. It is shown that this transformed GI/D/1
system has special structures which make the computation of the matrix R simple and efficient, thereby reducing the number of multiplications in each iteration significantly. As a result we were able to keep the computation time very low. Moreover, use of the resulting structural properties makes the computation of the distribution of queue length of the transformed system efficient. The computation of the distribution of waiting time is also shown to be simple by exploiting the special structures.


---

## Body

## 1. Introduction

In most communication systems studied in discrete time, the interarrival times of packets usually follow independent general distribution. Hence, this arrival of packets in an ATM switch can be represented as independent general distribution. On the other hand, multiple packets (53 byte cells) are transmitted simultaneously through identical transmission lines in an ATM switch. Therefore, an ATM switch can be considered as a multiserver system with deterministic service time as the packets are of same size. Hence, the performance of an ATM switch can be analyzed by studying a GI/D/k system with finite support on interarrival times.To the best of our knowledge GI/D/k systems in discrete time have never been analyzed with focus on computational efficiency in the literature. In this paper, matrix-geometric method is used to analyze a class of GI/D/k systems in discrete time with focus on computational efficiency. The class of GI/D/k systems considered in this paper has finite support for the interarrival times and the service is first-in first-out (FIFO) order. The idea used here to analyze the waiting time distribution of this multiserver system is due to Crommelin [1]. This idea is also used in analyzing the distribution of waiting time of multiserver systems with constant service time duration for both continuous and discrete times such as [2–20]. The idea is that the assignment of customers to servers in cyclic order does not change the distribution of waiting time in a multiserver system where customers are served in FIFO order in the identical servers with same constant service time. In this case, only an arbitrary server can be studied to compute the waiting time distribution. Crommelin [1] analyzed the waiting time distribution of an M/D/c system by converting the system to an Ec/D/1 system. The work of Crommelin [1] is followed by some studies on deterministic service time multiserver systems in both continuous time and discrete time. Iversen [15] decomposed an M/D/rk queue with FIFO into Ek/D/r queues with FIFO for analyzing waiting time distribution. Franx [8] developed expression for the waiting time distribution of the M/D/c queue by a full probabilistic analysis requiring neither generating functions nor Laplace transforms. Wittevrongel and Bruneel [19] used transform-based method to analyze multiserver system with constant service time and correlated arrival process. Nishimura [17] studied a MAP/D/N system in continuous time using spectral technique. Takine [18] analyzed a MAP/D/s system and computed the distribution of queue length by first characterizing the sojourn time distribution. Kim and Chaudhry [21] also computed the distribution of queue length of a discrete time GI/D/c queue after computing the distribution of waiting time by using distributional Little's law and transform-based method. Exact elapsed time from the time of last arrival is approximated with time average in [21]. Chaudhry et al. [7] analyzed the distribution of waiting time of a MAP/D/k system in discrete time. Later Alfa [3] gave a simple and more efficient computational scheme for the MAP/D/k system in discrete time. Alfa [2] also carried out algorithmic analysis for a discrete time BMAP/D/k system by exploiting its structural properties. Chaudhry et al. [7] and Alfa [2, 3] used matrix geometric method [22] for analysis. There is some work on the GI/D/c queues such as Wuyts and Bruneel [4] and Gao et al. [10–14] which used transform-based method. There are other algorithms for multiserver queues with constant service times such as GI/D/1 and GI/D/c queues by Chaudhry [23], GIX/D/c queue by Chaudhry and Kim [5], MX/D/c queue by Chaudhry et al. [6] and Franx [9], Ph/D/c queues by van Hoorn [16], and an Ek/D/r queue by Xerocostas and Demertzes [20].In this paper, the cyclic assignment of customers to servers is used to model a class of GI/D/k systems in discrete time as a single server system with the assumption of first-in first-out (FIFO) service order. The modeling as single server system sets up the distribution of queue length as a quasi-birth-death (QBD) which has some structural properties. Analysis of the GI/D/k system is carried out efficiently by exploiting these structural properties in this paper. This paper has three contributions—reductions in the computational complexities of (i) the matrix R, (ii) the distribution of queue length of the transformed system, and (iii) the distribution of waiting time. The first contribution is that the time complexity of the computation of the matrix R is reduced by decreasing the number of multiplications in each iteration from O(n3k3d3) to O(n3kd2) while requiring the same number of iteration as the natural iteration. Here, n is the support of general interarrival times, k is the number of servers, and d is the duration of service. The second contribution is that the distribution of queue length of the transformed system is computed efficiently by exploiting the special structures of the system. The third and most important contribution is that the computation of the distribution of waiting time is simplified. Chaudhry et al. [7] and Alfa [2, 3] computed the distribution of waiting time using the technique for phase type service time distribution. However, it is shown in this paper that the computation of the distribution of waiting time for deterministic service time distribution in discrete time does not need the complicated steps of phase type service time distribution.The rest of the paper is organized as follows. Section2 introduces the GI/D/k system and explains the modeling of this multiserver system into a single server system. Section 3 discusses the special structures of the matrices of our model and exploitation of those structures for efficient computation of the matrix R. The computation of the distributions of queue length and waiting time is explained in Sections 4 and 5, respectively. Some numerical examples are provided in Section 6 to show the suitability of our proposed method to compute the matrix R. Section 7 concludes the paper.


## 2. The GI/D/k System

The class of GI/D/k system in discrete time considered here has finite support for interarrival time which is of general distribution. Moreover, the customers are served in FIFO service order in the identical servers which have constant service time.The interarrival times,A, of this multiserver system are general, independent, and identically distributed (iid) with distribution av=Pr{A=v}, v=1,2,….We letλ=1/E[A] be the mean arrival rate of customers to the system with 0<λ≤1. The interarrival times can be represented as a Markov chain where each state represents the elapsed time using the approach followed in the analysis of discrete time GI/G/1 and GIX/G/1 systems by Alfa [24]. General independent interarrival times can be represented as an absorbing Markov chain with a transition matrix PA=[Tt01,n1] where T is a square matrix of order n, t is a column vector of order n, t=1n-T1n, 0i,j is an i×j matrix of zeros, and 1i is a column vector of ones of order i. The case of infinite interarrival times can be captured for n=∞. However, the interarrival times are finite if data is collected from practical systems. Therefore, general independent interarrival times with finite support (n<∞) are considered in this paper. The probabilities of interarrival times can be represented by a vector a̲=[a1,a2,…,an]. Here, a̲1n=1 and a0=0. This general arrival process can be defined by phase type distribution (α,T) of dimension n in elapsed time representation. The parameters α, T, and t are given as α=[1,0,0…,0], Ti,j=(1-∑j=0iaj)/(1-∑v=0i-1av)=ãi, for 1≤i<n, j=i+1, Ti,j=0, for all j≠i+1 and t=[ã̃1ã̃2⋯ã̃n-11]′ where ã̃i=1-ãi for 1≤i≤n-1 and W′ is the transpose of matrix W.The general independent arrival process can be represented by two square matricesD0 and D1 of order n where (Dl)ij(l=0,1and1≤i,j≤n) represents a transition from phase i to phase j with l arrivals. Now, D0 and D1 can be represented in terms of phase type distribution (α,T). Here, D0=T and D1=tα=[t0n,n-1]. The first column of D1 is nonzero and this column is the vector t. The matrices D0 and D1 are similar to the matrices of zero or one arrival for MAP distributed interarrival times. However, they are special cases because the matrices D0 and D1 only capture general and independently distributed interarrival times.A discrete time Markov chain{(N(t),J(t));t=0,1,2,…;N(t)∈{0,1,2,…};J(t)∈{1,2,3,…,n}} can be considered to represent this general arrival process where N(t) is the number of arrivals up to and including time t and J(t) is the phase of the next arrival at time t.The matrixD=D0+D1 is stochastic. If π̃=[π̃1,π̃2,…,π̃n] is the invariant probability vector for arrival in n different phases then π̃=π̃D and π̃1n=1. The arrival rate is given by λ=π̃D11n=π̃t.There arek identical servers in a GI/D/k system. The service times are of constant duration d units (d≥1). Service time can be represented as phase type distribution in elapsed time form with representation (β,S) of dimension d where β=[1,0,0,…,0] and S=[0d-1,1Id-1001,d-1], Ij is an identity matrix of order j. Another column vector s=1d-S1d can be defined for this phase type distribution (β,S).The condition that ensures the stability of GI/D/k system is ρ<1 where ρ=λd/k. We assume that the system is stable (i.e., ρ<1).


### 2.1. The Model of Arrivals to an Arbitrary Server

Let us assume that customers are served in first-in first-out (FIFO) order in this GI/D/k system. It can also be assumed that the servers are assigned with customers in a cyclic manner particularly during the idle times. This assignment of customers to servers does not affect the waiting times experienced by the customers as the servers are identical. Using the same approach of Alfa [2, 3] we can consider the case of the jth (1≤j≤k) server. In this approach if the first customer is assigned to the first server then the jth server will be assigned with j, j+k, j+2k,…th customers. In this case we only need to study one server, an arbitrary server j  (1≤j≤k). This arbitrary server sees the arrival of customers to be according to another independent general process which we call GIk. This GIk arrival process is a k-fold convolution of original arrival process to the multiserver system. GIk can be defined by two square matrices C0 and C1 of order nk. The element (Cv)ij (v=0,1 and 1≤i,j≤nk) represents transition from phase i to phase j with v arrivals. Both C0 and C1 can be considered to comprise of k2 square block matrices of order n. In this representation (C0)i,j=D0 for i=j,1≤i≤k, (C0)i,j=D1 for j=i+1,1≤i<k and (C0)i,j=0n,n for 1≤i≤k and i≠j or j≠i+1. Similarly (C1)i,j=D1 for i=k,j=1 and (C1)i,j=0n,n in all other cases for 1≤i,j≤k. The first column and the last row of C0 are zero as the first column and last row of D0 are zero. Only the last n rows of C1 are nonzero and only the first element of each of these n rows is nonzero.The matricesC0 and C1 exhibit behavior similar to the matrices D0 and D1 except that Cv (v=0,1) is a new matrix representing the arrival of v supercustomers where a supercustomer is the last customer to form a group of k customers or the arrival of v regular customers to an arbitrary server. On the other hand, Dv  (v=0,1) represents arrival of v regular customers to the GI/D/k system. Now, this new arrival process GIk to an arbitrary server can be represented by a discrete time Markov chain {(N̂(t),Ĵ(t));t=0,1,2,…;N̂(t)∈{0,1,2,…};Ĵ(t)∈{1,2,…,nk}} where N̂(t) is the number of supercustomers that have arrived by time t and Ĵ(t) is the phase of the next arrival at time t.The matrixC=C0+C1 is stochastic. If δ, a 1×nk matrix, is the invariant probability vector for arrival in different nk phases then δ=δC and δ1nk=1. δ can also be represented as δ=[δ1,δ2,δ3,…,δk] where δi(1≤i≤k) is of size 1×n. The supercustomer arrival rate is given by λ*=λ/k=δC11nk=δkt.


### 2.2. The New ModelGkI/D/1 System

A model can be developed for a single server queue with arrivalGIk and deterministic service time of d units following the techniques of Alfa [3]. This single server queue is called GIk/D/1 system and the resulting model is of QBD process type for which standard methods can be used [22]. Let us assume that Lt, Jt, and St represent the number of supercustomers, the arrival phase of the next supercustomer, and the elapsed time of the supercustomer in service, respectively, at time t in this QBD model. The state space can be represented by {Lt,Jt,St} for Lt≥1 and the state space is {0,Jt} for Lt=0. The total state space is {(0,j)∪(i,j,l):i≥1,1≤j≤nk,1≤l≤d}. The transition matrix for this GIk/D/1 system is(2.1)P=[B0B1EA1A0A2A1A0A2A1A0⋱⋱⋱],
where B0=C0, B1=C1⊗β, E=C0⊗s, A2=C0⊗(sβ), A1=C0⊗S+C1⊗(sβ), and A0=C1⊗S.Letx=[x0,x1,x2,…] be the stationary distribution of P where x=xP, x01nk+∑i=1∞xi1nkd=1 and xj (j≥0) is the probability vector representing j supercustomers in the GIk/D/1 system. During the computation of the stationary distribution of queue length x, a matrix R is computed first which is followed by the computations of x0 and x1. Then xi(i≥2) is computed using the matrix-geometric result xi=xi-1R=x1Ri-1. Here, the matrix R is of order nkd×nkd which is the minimal nonnegative solution of R=A0+RA1+R2A2. The entry Ru,v(1≤u,v≤nkd) of R represents the expected number of visits into (m+1,v), starting from (m,u), before the first return to level m (m≥1). Here, m is the number of supercustomers in the system and u and v are any of the nkd phases of the system.


## 2.1. The Model of Arrivals to an Arbitrary Server

Let us assume that customers are served in first-in first-out (FIFO) order in this GI/D/k system. It can also be assumed that the servers are assigned with customers in a cyclic manner particularly during the idle times. This assignment of customers to servers does not affect the waiting times experienced by the customers as the servers are identical. Using the same approach of Alfa [2, 3] we can consider the case of the jth (1≤j≤k) server. In this approach if the first customer is assigned to the first server then the jth server will be assigned with j, j+k, j+2k,…th customers. In this case we only need to study one server, an arbitrary server j  (1≤j≤k). This arbitrary server sees the arrival of customers to be according to another independent general process which we call GIk. This GIk arrival process is a k-fold convolution of original arrival process to the multiserver system. GIk can be defined by two square matrices C0 and C1 of order nk. The element (Cv)ij (v=0,1 and 1≤i,j≤nk) represents transition from phase i to phase j with v arrivals. Both C0 and C1 can be considered to comprise of k2 square block matrices of order n. In this representation (C0)i,j=D0 for i=j,1≤i≤k, (C0)i,j=D1 for j=i+1,1≤i<k and (C0)i,j=0n,n for 1≤i≤k and i≠j or j≠i+1. Similarly (C1)i,j=D1 for i=k,j=1 and (C1)i,j=0n,n in all other cases for 1≤i,j≤k. The first column and the last row of C0 are zero as the first column and last row of D0 are zero. Only the last n rows of C1 are nonzero and only the first element of each of these n rows is nonzero.The matricesC0 and C1 exhibit behavior similar to the matrices D0 and D1 except that Cv (v=0,1) is a new matrix representing the arrival of v supercustomers where a supercustomer is the last customer to form a group of k customers or the arrival of v regular customers to an arbitrary server. On the other hand, Dv  (v=0,1) represents arrival of v regular customers to the GI/D/k system. Now, this new arrival process GIk to an arbitrary server can be represented by a discrete time Markov chain {(N̂(t),Ĵ(t));t=0,1,2,…;N̂(t)∈{0,1,2,…};Ĵ(t)∈{1,2,…,nk}} where N̂(t) is the number of supercustomers that have arrived by time t and Ĵ(t) is the phase of the next arrival at time t.The matrixC=C0+C1 is stochastic. If δ, a 1×nk matrix, is the invariant probability vector for arrival in different nk phases then δ=δC and δ1nk=1. δ can also be represented as δ=[δ1,δ2,δ3,…,δk] where δi(1≤i≤k) is of size 1×n. The supercustomer arrival rate is given by λ*=λ/k=δC11nk=δkt.


## 2.2. The New ModelGkI/D/1 System

A model can be developed for a single server queue with arrivalGIk and deterministic service time of d units following the techniques of Alfa [3]. This single server queue is called GIk/D/1 system and the resulting model is of QBD process type for which standard methods can be used [22]. Let us assume that Lt, Jt, and St represent the number of supercustomers, the arrival phase of the next supercustomer, and the elapsed time of the supercustomer in service, respectively, at time t in this QBD model. The state space can be represented by {Lt,Jt,St} for Lt≥1 and the state space is {0,Jt} for Lt=0. The total state space is {(0,j)∪(i,j,l):i≥1,1≤j≤nk,1≤l≤d}. The transition matrix for this GIk/D/1 system is(2.1)P=[B0B1EA1A0A2A1A0A2A1A0⋱⋱⋱],
where B0=C0, B1=C1⊗β, E=C0⊗s, A2=C0⊗(sβ), A1=C0⊗S+C1⊗(sβ), and A0=C1⊗S.Letx=[x0,x1,x2,…] be the stationary distribution of P where x=xP, x01nk+∑i=1∞xi1nkd=1 and xj (j≥0) is the probability vector representing j supercustomers in the GIk/D/1 system. During the computation of the stationary distribution of queue length x, a matrix R is computed first which is followed by the computations of x0 and x1. Then xi(i≥2) is computed using the matrix-geometric result xi=xi-1R=x1Ri-1. Here, the matrix R is of order nkd×nkd which is the minimal nonnegative solution of R=A0+RA1+R2A2. The entry Ru,v(1≤u,v≤nkd) of R represents the expected number of visits into (m+1,v), starting from (m,u), before the first return to level m (m≥1). Here, m is the number of supercustomers in the system and u and v are any of the nkd phases of the system.


## 3. Structures of the Matrices

This section explains the structures of the block-matricesB0, B1, E, A2, A1, and A0 and how their structures are exploited for efficient computation of the matrix R. The structures of these matrices are also fully exploited in computing the stationary distribution of the queue length in Section 4.Here,B0 and B1 are matrices of order nk×nk and nk×nkd, respectively. The first column and the last row of the matrix B0 are zero. On the other hand, only the last n rows of the matrix B1 are nonzero and among these last n rows only the first column is nonzero which is t. The matrix B1 can be represented as B1=[0n(k-1),nd0n(k-1),nd(k-1)L̂0n,nd(k-1)] where L̂=[t0n,nd-1].The matrixE is of order nkd×nk which can be represented as
(3.1)E=[L̃M̃L̃M̃⋱⋱L̃M̃L̃],
where
(3.2)L̃=[ã1sã2s⋱ãn-1s],M̃=[ã̃1sã̃2s⋮ã̃n-1ss].L̃=D0⊗s is a matrix of order nd×n. Moreover, L̃ has (n-1) nonzero columns and each column contains only one nonzero element. L̃ has (n-1) nonzero elements in (jd,j+1)th (1≤j≤n-1) positions. M̃=D1⊗s is also a matrix of order nd×n matrix. Furthermore, only the first column of M̃ is nonzero which contains n nonzero elements. M̃ has nonzero elements in (jd,1)th (1≤j≤n) positions.The matricesA2, A1, and A0 are of size nkd×nkd. The matrix A2 can be represented as(3.3)A2=[L̃1M̃1L̃1M̃1⋱⋱L̃1M̃1L̃1],
where L̃1 and M̃1 are both matrices of order nd×nd, L̃1=D0⊗(sβ) and M̃1=D1⊗(sβ). Hence,(3.4)L̃1=[ã1(sβ)ã2(sβ)⋱ãn-1(sβ)],M̃1=[ã̃1(sβ)ã̃2(sβ)⋮ã̃n-1(sβ)(sβ)].There are(n-1) nonzero columns in L̃1 and each column contains only one nonzero element. L̃1 has these nonzero elements in (jd,jd+1)th (1≤j≤n-1) positions. Only the first column of M̃1 is nonzero with n nonzero elements in (jd,1)th (1≤j≤n) positions.The matrixA1 can be represented as
(3.5)A1=[L̃2M̃2L̃2M̃2⋱⋱L̃2M̃2M̃1L̃2],
where
(3.6)L̃2=[ã1Sã2S⋱ãn-1S],M̃2=[ã̃1Sã̃2S⋮ã̃n-1SS].Here,L̃2 and M̃2 are both matrices of order nd×nd. L̃2=D0⊗S and M̃2=D1⊗S. There are (n-1)(d-1) nonzero columns, (n-1)(d-1) nonzero rows, and (n-1)(d-1) nonzero elements in L̃2. Moreover, each nonzero row or each nonzero column of L̃2 contains only one nonzero element. On the other hand, there are (d-1) nonzero columns in M̃2 from column 2 to d each having n nonzero elements.The matrixA0 can be represented as A0=[M̃2]. A0 has only n(d-1) nonzero rows which are in its last nd rows and (d-1) nonzero columns from column 2 to column d.


### 3.1. Computation of the MatrixR

The matrixR can be computed using any of the three linearly convergent formulae from Neuts [22]—R(i+1):=A0+R(i)A1+R2(i)A2, R(i+1):=(A0+R2(i)A2)(Inkd-A1)-1, and R(i+1):=A0(Inkd-A1-R(i)A2)-1 which are known as natural, traditional, and U-based iterations, respectively. Here, R(i) is the ith iteration value of R, R(0):=0nkd,nkd, U is a square matrix of the same order as the matrix R and U=A1+RA2. The matrix U is the minimal nonnegative solution of U=A1+A0(I-U)-1A2. The entry Uu,v(1≤u,v≤nkd) represents the taboo probability of starting from (m,u) for m≥1 the eventual visit of the Markov chain to level m by visiting (m,v) under the taboo of level (m-1). Any of these three methods terminate in the ith (i≥1) iteration for ∥R(i)-R(i-1)∥∞<ε where ε is a very small positive constant. In each step, natural itearation requires two matrix multiplications whereas traditional iteration requires three matrix multiplications and one matrix inversion. On the other hand, U-based iteration requires two matrix multiplications and one matrix inversion in each step. Traditional iteration requires less number of iterations than natural iteration does while the number of iterations required for U-based iteration is less than the numbers of iterations for other two methods. Alfa and Xue [25] developed efficient inversion technique for U-based iteration for the matrix R in analyzing the queue length of a GI/G/1 system in discrete time.There are several efficient methods for computing the matrixR for a QBD such as logarithmic reduction algorithm of Latouche and Ramaswami [26], cyclic reduction technique of Bini and Meini [27], and invariant subspace approach of Akar and Sohraby [28]. These quadratically convergent algorithms compute a matrix G which is the minimal nonnegative solution to G=A2+A1G+A0G2. Here, G is an m×m matrix for m×m matrices A2, A1, and A0 where Gi,j  (1≤i,j≤m) represents the probability that starting from state (N+1,i) the Markov chain eventually visits the level N  (N≥1) and does so by visiting the state (N,j). The matrices R, U, and G are related and these relations can be found in [26]. Similar to the computation of the matrix R, the matrix G can be computed using three different iterations—natural, traditional, and U-based. Generally a zero matrix is used for the initial value for iterations of the matrix G (i.e., G=0m,m). Latouche and Ramaswami [26] showed that if SG(ε) iterative steps are needed to converge for a very small positive constant ε such that ∥G(SG(ε))-G(SG(ε)-1)∥∞<ε for the (l+1)-st natural iteration G(l+1)=A2+A1G(l)+A0G2(l), then logarithmic reduction requires not more than log2SG(ε) iterations and overall complexity of their algorithm is O(m3log2SG(ε)). They also developed logarithmic reduction technique for the matrix R. However, the quadratic convergent algorithms [26–28] are not considered here for computing the matrix R as these algorithms involve inversion of matrices which cannot exploit the special structures of the matrices A2, A1, and A0.The complexity of computing the matrixR can be reduced by using the structure of A0. The first nd(k-1) rows are zero in A0. Therefore, the first nd(k-1) rows are zero in R. Among the last nd rows of A0, n rows are zero which are rows nd(k-1)+d, nd(k-1)+2d, nd(k-1)+3d,…,nkd. The corresponding rows are zero for R. The matrix R can be represented as follows:(3.7)R=[R1R2R3⋯Rk].
Here, Rj(1≤j≤k) is a matrix of order nd×nd which contains d,2d,3d,…,ndth rows as zero rows.The computations ofR(i)A2 and R(i)A1 require (2k-1) and 2k block-matrix multiplications, respectively, where each block matrix is of order nd×nd. Similarly, the computation of R2(i)A2 requires multiplication of R(i) and R(i)A2 which involves k block-matrix multiplications. Inversions of both the matrices (Inkd-A1) and (Inkd-A1-R(i)A2) generate k2 block matrices each of order nd×nd. Inversion of the matrix (Inkd-A1) needs to be computed only once for traditional iteration whereas inversion of the matrix (Inkd-A1-R(i)A2) needs to be computed in each iteration for U-based iteration. On the other hand, the multiplication of (A0+R2(i)A2) and (Inkd-A1)-1 in traditional iteration and the multiplication of A0 and (Inkd-A1-R(i)A2)-1 in U-based iteration involve k2 and k block-matrix multiplications, respectively, where each block matrix is of order nd×nd. Each iteration of natural method requires substantially less multiplications than each step of traditional and U-based methods do. Therefore, we use natural iteration for developing our algorithm for computing the matrix R in this paper.The matrixRj(1≤j≤k) of (3.7) can be computed using the natural iteration with the following equation:(3.8)Rj(i+1):={M̃2+R1(i)L̃2+Rk(i)M̃1+Rk(i)R1(i)L̃1,j=1,Rj-1(i)M̃2+Rj(i)L̃2+Rk(i)Rj-1(i)M̃1+Rk(i)Rj(i)L̃1,1<j≤k,
where Rj(i)(i≥0,1≤j≤k) is the ith iteration value of Rj with the initial value Rj(0):=0nd,nd(1≤j≤k). The terminating condition for this iteration is same as the terminating condition for natural, traditional, and U-based methods but needs to consider only n(d-1) nonzero rows of (R(i)-R(i-1)).The number of multiplications required to computeRj(i)L̃2(1≤j≤k) is n(n-1)(d-1)2. Similarly, the number of multiplications required to compute both Rj(i)L̃1 and Rj(i)M̃1 for 1≤j≤k is n(n-1)(d-1). In the same way, the number of multiplications required to compute Rj-1(i)M̃2(1<j≤k) is n(n-1)(d-1)2. On the other hand, the numbers of multiplications required to compute Rk(i)Rj-1(i)M̃1(1<j≤k) by block-matrix multiplication of Rk(i) and Rj-1(i)M̃1(1<j≤k) and Rk(i)Rj(i)L̃1(1≤j≤k) by block-matrix multiplication of Rk(i) and Rj(i)L̃1(1≤j≤k) are n2(d-1)2 and n2(n-1)(d-1)2, respectively. Hence, the numbers of multiplications required for computing R1(i+1) from R1(i) and Rj(i+1)(2≤j≤k) from Rj(i) in one iteration are n(n-1)(d-1)(d+1)+n2(n-1)(d-1)2 and n3(d-1)2+2nd(n-1)(d-1), respectively. Therefore, the total number of multiplications required to compute all k block matrices Rj(i+1)'s (1≤j≤k) of the matrix R(i+1) in one iteration is n2(n-1)(d-1)2+n(n-1)(d-1)(d+1)+(k-1)(n3(d-1)2+2nd(n-1)(d-1))=O(n3kd2). By exploiting the special structures of the matrices the computational complexity is decreased in each iteration from O(n3k3d3) to O(n3kd2) which is substantial reduction in computations. Moreover, the memory requirement for our method is also less than that of natural, traditional, U-based methods and quadratic convergent algorithms.


## 3.1. Computation of the MatrixR

The matrixR can be computed using any of the three linearly convergent formulae from Neuts [22]—R(i+1):=A0+R(i)A1+R2(i)A2, R(i+1):=(A0+R2(i)A2)(Inkd-A1)-1, and R(i+1):=A0(Inkd-A1-R(i)A2)-1 which are known as natural, traditional, and U-based iterations, respectively. Here, R(i) is the ith iteration value of R, R(0):=0nkd,nkd, U is a square matrix of the same order as the matrix R and U=A1+RA2. The matrix U is the minimal nonnegative solution of U=A1+A0(I-U)-1A2. The entry Uu,v(1≤u,v≤nkd) represents the taboo probability of starting from (m,u) for m≥1 the eventual visit of the Markov chain to level m by visiting (m,v) under the taboo of level (m-1). Any of these three methods terminate in the ith (i≥1) iteration for ∥R(i)-R(i-1)∥∞<ε where ε is a very small positive constant. In each step, natural itearation requires two matrix multiplications whereas traditional iteration requires three matrix multiplications and one matrix inversion. On the other hand, U-based iteration requires two matrix multiplications and one matrix inversion in each step. Traditional iteration requires less number of iterations than natural iteration does while the number of iterations required for U-based iteration is less than the numbers of iterations for other two methods. Alfa and Xue [25] developed efficient inversion technique for U-based iteration for the matrix R in analyzing the queue length of a GI/G/1 system in discrete time.There are several efficient methods for computing the matrixR for a QBD such as logarithmic reduction algorithm of Latouche and Ramaswami [26], cyclic reduction technique of Bini and Meini [27], and invariant subspace approach of Akar and Sohraby [28]. These quadratically convergent algorithms compute a matrix G which is the minimal nonnegative solution to G=A2+A1G+A0G2. Here, G is an m×m matrix for m×m matrices A2, A1, and A0 where Gi,j  (1≤i,j≤m) represents the probability that starting from state (N+1,i) the Markov chain eventually visits the level N  (N≥1) and does so by visiting the state (N,j). The matrices R, U, and G are related and these relations can be found in [26]. Similar to the computation of the matrix R, the matrix G can be computed using three different iterations—natural, traditional, and U-based. Generally a zero matrix is used for the initial value for iterations of the matrix G (i.e., G=0m,m). Latouche and Ramaswami [26] showed that if SG(ε) iterative steps are needed to converge for a very small positive constant ε such that ∥G(SG(ε))-G(SG(ε)-1)∥∞<ε for the (l+1)-st natural iteration G(l+1)=A2+A1G(l)+A0G2(l), then logarithmic reduction requires not more than log2SG(ε) iterations and overall complexity of their algorithm is O(m3log2SG(ε)). They also developed logarithmic reduction technique for the matrix R. However, the quadratic convergent algorithms [26–28] are not considered here for computing the matrix R as these algorithms involve inversion of matrices which cannot exploit the special structures of the matrices A2, A1, and A0.The complexity of computing the matrixR can be reduced by using the structure of A0. The first nd(k-1) rows are zero in A0. Therefore, the first nd(k-1) rows are zero in R. Among the last nd rows of A0, n rows are zero which are rows nd(k-1)+d, nd(k-1)+2d, nd(k-1)+3d,…,nkd. The corresponding rows are zero for R. The matrix R can be represented as follows:(3.7)R=[R1R2R3⋯Rk].
Here, Rj(1≤j≤k) is a matrix of order nd×nd which contains d,2d,3d,…,ndth rows as zero rows.The computations ofR(i)A2 and R(i)A1 require (2k-1) and 2k block-matrix multiplications, respectively, where each block matrix is of order nd×nd. Similarly, the computation of R2(i)A2 requires multiplication of R(i) and R(i)A2 which involves k block-matrix multiplications. Inversions of both the matrices (Inkd-A1) and (Inkd-A1-R(i)A2) generate k2 block matrices each of order nd×nd. Inversion of the matrix (Inkd-A1) needs to be computed only once for traditional iteration whereas inversion of the matrix (Inkd-A1-R(i)A2) needs to be computed in each iteration for U-based iteration. On the other hand, the multiplication of (A0+R2(i)A2) and (Inkd-A1)-1 in traditional iteration and the multiplication of A0 and (Inkd-A1-R(i)A2)-1 in U-based iteration involve k2 and k block-matrix multiplications, respectively, where each block matrix is of order nd×nd. Each iteration of natural method requires substantially less multiplications than each step of traditional and U-based methods do. Therefore, we use natural iteration for developing our algorithm for computing the matrix R in this paper.The matrixRj(1≤j≤k) of (3.7) can be computed using the natural iteration with the following equation:(3.8)Rj(i+1):={M̃2+R1(i)L̃2+Rk(i)M̃1+Rk(i)R1(i)L̃1,j=1,Rj-1(i)M̃2+Rj(i)L̃2+Rk(i)Rj-1(i)M̃1+Rk(i)Rj(i)L̃1,1<j≤k,
where Rj(i)(i≥0,1≤j≤k) is the ith iteration value of Rj with the initial value Rj(0):=0nd,nd(1≤j≤k). The terminating condition for this iteration is same as the terminating condition for natural, traditional, and U-based methods but needs to consider only n(d-1) nonzero rows of (R(i)-R(i-1)).The number of multiplications required to computeRj(i)L̃2(1≤j≤k) is n(n-1)(d-1)2. Similarly, the number of multiplications required to compute both Rj(i)L̃1 and Rj(i)M̃1 for 1≤j≤k is n(n-1)(d-1). In the same way, the number of multiplications required to compute Rj-1(i)M̃2(1<j≤k) is n(n-1)(d-1)2. On the other hand, the numbers of multiplications required to compute Rk(i)Rj-1(i)M̃1(1<j≤k) by block-matrix multiplication of Rk(i) and Rj-1(i)M̃1(1<j≤k) and Rk(i)Rj(i)L̃1(1≤j≤k) by block-matrix multiplication of Rk(i) and Rj(i)L̃1(1≤j≤k) are n2(d-1)2 and n2(n-1)(d-1)2, respectively. Hence, the numbers of multiplications required for computing R1(i+1) from R1(i) and Rj(i+1)(2≤j≤k) from Rj(i) in one iteration are n(n-1)(d-1)(d+1)+n2(n-1)(d-1)2 and n3(d-1)2+2nd(n-1)(d-1), respectively. Therefore, the total number of multiplications required to compute all k block matrices Rj(i+1)'s (1≤j≤k) of the matrix R(i+1) in one iteration is n2(n-1)(d-1)2+n(n-1)(d-1)(d+1)+(k-1)(n3(d-1)2+2nd(n-1)(d-1))=O(n3kd2). By exploiting the special structures of the matrices the computational complexity is decreased in each iteration from O(n3k3d3) to O(n3kd2) which is substantial reduction in computations. Moreover, the memory requirement for our method is also less than that of natural, traditional, U-based methods and quadratic convergent algorithms.


## 4. Stationary Distribution of P forGkI/D/1 System

Two methods for computing stationary distributionx for the Markov chain represented by the matrix P are briefly presented here which are followed by Section 4.1. Section 4.1 explains how the second method can be made efficient by utilizing the special structures of the matrices.In the first method an invariant probability vectorz=[z0,z1] needs to be computed for a stochastic matrix P̃=[B0B1EA1+RA2]. The row vector z is computed from zP̃=z and z01nk+z11nkd=1. Then x0 and x1 can be expressed as x0=v-1z0 and x1=v-1z1 where v=z01nk+z1(Inkd-R)-11nkd.In the second method, we need to computeh̃, the invariant probability vector of a stochastic matrix H defined as H=A1+RA2+E(Ink-B0)-1B1. Next a positive real number h is computed as follows h=h̃(E(Ink-B0)-11nk+(Inkd-R)-11nkd). On the other hand, h̃ is computed using the relations h̃=h̃H and h̃1nkd=1. x0 and x1 can be computed as x1=h-1h̃, x0=x1E(Ink-B0)-1.


### 4.1. Exploiting Structures for Stationary Distribution

The stationary distribution ofx can be computed by first computing the probability vector x1 which can be made efficient by exploiting the special structures of the matrices H, E(Ink-B0)-11nk and (Inkd-R)-11nkd. The special structures of the matrices (Inkd-R)-1, (Inkd-R)-11nkd, (Ink-B0)-1, E(Ink-B0)-1, E(Ink-B0)-11nk, E(Ink-B0)-1B1, and H are due to the special structures of the matrices B0, B1, E, and R. Here, how the computations of these matrices can be carried out efficiently is explained.Computing(Inkd-R)-1=[Ind(k-1)0nd(k-1),ndΦΥ] requires O(n3kd3) multiplications by using block-matrix inversion [29] where Υ=R̃k-1, Φ=ΥΨ, Ψ=[R1,R2,R3,…,Rk-1], and R̃k=Ind-Rk. Φ can be expressed as Φ=[Φ1,Φ2,Φ3,…,Φk-1] where Φj=ΥRj(1≤j≤k-1).(Ink-B0)-1 can be represented as(4.1)(Ink-B0)-1=[V1V2⋯Vk-1VkV1⋯Vk-2Vk-1⋱⋮⋮V1V2V1],
where V1=D̃0-1, Vi=(D̃0-1D1)i-1D̃0-1(1<i≤k), and D̃0=In-D0. The computation of (Ink-B0)-1 requires O(n2k+nk2) instead of O(n3k3) multiplications by utilizing the special structure of (Ink-B0).The matricesE(Ink-B0)-1 and E(Ink-B0)-1B1 can be represented as(4.2)E(Ink-B0)-1=[Ṽ1Ṽ2…Ṽk-1ṼkṼ1…Ṽk-2Ṽk-1⋱⋮⋮Ṽ1Ṽ2Ṽ1],E(Ink-B0)-1B1=[Ṽ̃1Ṽ̃2⋮Ṽ̃k-1Ṽ̃k],
where Ṽ1=L̃V1, Ṽi=L̃Vi+M̃Vi-1, for 1<i≤k, and Ṽ̃i=Ṽk-i+1L̂(1≤i≤k). Here, Ṽi and Ṽ̃i for 1≤i≤k are matrices of order nd×n and nd×nd, respectively. Here, multiplying E by (Ink-B0)-1 and E(Ink-B0)-1 by B1 both require O(n2k) multiplications due to utilizing special structures of matrices. If there were no special structures of E, (Ink-B0)-1 and B1 then the multiplying E by (Ink-B0)-1 and E(Ink-B0)-1 by B1 would require O(n3k3d) and O(n3k2d3) multiplications, respectively.The matrixH of can be represented in the following form:(4.3)H=[H1,1H1,2H2,1H2,2H2,3H3,1H3,3H3,4⋮⋱⋱Hk-1,1Hk-1,k-1Hk-1,kHk,1Hk,2Hk,3⋯Hk,k-1Hk,k],
where Hi,j(1≤i,j≤k) is a matrix of order nd×nd. Here,(4.4)Hi,j={L̃2+Ṽ̃1,i=j=1,M̃2,1≤i≤k-1,j=i+1,L̃2,1<i=j<k,Ṽ̃i,1<i<k,j=1,R1L̃1+Ṽ̃k,i=k,j=1,Ri-1M̃1+RiL̃1,i=k,1<j<k,Rk-1M̃1+RkL̃1+L̃2,i=j=k.We can recursively defineUk=Ũk=Hk-1,k(I-Hk,k)-1, Ui=Hi-1,i(I-Hi,i-Ũi+1Hk,i)-1 for 2≤i≤k-1, Ũi=UiŨi+1 for 2≤i≤k-1, Ũ̃2=U2, Ũ̃i=Ũ̃i-1Ui for 3≤i≤k, and U1=H1,1+∑i=2kŨ̃iHi,1.Let us define another invariant probability vectorũ=[ũ1,ũ2,ũ3,…,ũk] of the stochastic matrix H where ũi(1≤i≤k) is a row vector having nd scalar elements. Here, ũ is computed from ũ1U1=ũ1, ũ11nd=1 and ũi=ũi-1Ui for 2≤i≤k. Next a positive real number u is computed as u=ũ(η+ζ) where η=E(Ink-B0)-11nk and ζ=(Inkd-R)-11nkd. η and ζ are both of size nkd×1. The column vectors η and ζ can be expressed as η=[η1,η2,η3,…,ηk-1,ηk] and ζ=[ζ1,ζ2,ζ3,…,ζk-1,ζk], respectively. Both ηi and ζi(1≤i≤k) are of size nd×1. ηk=Ṽ11n and ηj=ηj+1+Ṽk-j+11n for (1≤j≤k-1). There are only (n-1) nonzero elements in ηk and n nonzero elements in ηj for 1≤j≤k-1. On the other hand, ζi=1nd for 1≤i≤k-1 and ζk=(∑j=1k-1Φj+Υ)1nd.x1 can be computed as x1=u-1ũ using O(n3kd3) multiplications instead of O(n3k3d3) multiplications using the special structure of H. Let us partition x0 and x1 as x0=[κ1,κ2,κ3,…,κk] and x1=[ν1,ν2,ν3,…,νk], respectively, for efficient computation of x0. Here, κi and νi for 1≤i≤k are of size 1×n and 1×nd, respectively. κi(1≤i≤k) is computed as κi=∑j=1iνi-j+1Ṽj and these vector-matrix multiplications can be made efficient by considering only nonzero rows of Ṽj(1≤j≤k).xi(i≥2) can be computed as xi=xi-1R for (i≥2) using n(d-1)×nkd=n2kd(d-1) multiplications instead of n2k2d2 multiplications as R has only n(d-1) nonzero rows.


## 4.1. Exploiting Structures for Stationary Distribution

The stationary distribution ofx can be computed by first computing the probability vector x1 which can be made efficient by exploiting the special structures of the matrices H, E(Ink-B0)-11nk and (Inkd-R)-11nkd. The special structures of the matrices (Inkd-R)-1, (Inkd-R)-11nkd, (Ink-B0)-1, E(Ink-B0)-1, E(Ink-B0)-11nk, E(Ink-B0)-1B1, and H are due to the special structures of the matrices B0, B1, E, and R. Here, how the computations of these matrices can be carried out efficiently is explained.Computing(Inkd-R)-1=[Ind(k-1)0nd(k-1),ndΦΥ] requires O(n3kd3) multiplications by using block-matrix inversion [29] where Υ=R̃k-1, Φ=ΥΨ, Ψ=[R1,R2,R3,…,Rk-1], and R̃k=Ind-Rk. Φ can be expressed as Φ=[Φ1,Φ2,Φ3,…,Φk-1] where Φj=ΥRj(1≤j≤k-1).(Ink-B0)-1 can be represented as(4.1)(Ink-B0)-1=[V1V2⋯Vk-1VkV1⋯Vk-2Vk-1⋱⋮⋮V1V2V1],
where V1=D̃0-1, Vi=(D̃0-1D1)i-1D̃0-1(1<i≤k), and D̃0=In-D0. The computation of (Ink-B0)-1 requires O(n2k+nk2) instead of O(n3k3) multiplications by utilizing the special structure of (Ink-B0).The matricesE(Ink-B0)-1 and E(Ink-B0)-1B1 can be represented as(4.2)E(Ink-B0)-1=[Ṽ1Ṽ2…Ṽk-1ṼkṼ1…Ṽk-2Ṽk-1⋱⋮⋮Ṽ1Ṽ2Ṽ1],E(Ink-B0)-1B1=[Ṽ̃1Ṽ̃2⋮Ṽ̃k-1Ṽ̃k],
where Ṽ1=L̃V1, Ṽi=L̃Vi+M̃Vi-1, for 1<i≤k, and Ṽ̃i=Ṽk-i+1L̂(1≤i≤k). Here, Ṽi and Ṽ̃i for 1≤i≤k are matrices of order nd×n and nd×nd, respectively. Here, multiplying E by (Ink-B0)-1 and E(Ink-B0)-1 by B1 both require O(n2k) multiplications due to utilizing special structures of matrices. If there were no special structures of E, (Ink-B0)-1 and B1 then the multiplying E by (Ink-B0)-1 and E(Ink-B0)-1 by B1 would require O(n3k3d) and O(n3k2d3) multiplications, respectively.The matrixH of can be represented in the following form:(4.3)H=[H1,1H1,2H2,1H2,2H2,3H3,1H3,3H3,4⋮⋱⋱Hk-1,1Hk-1,k-1Hk-1,kHk,1Hk,2Hk,3⋯Hk,k-1Hk,k],
where Hi,j(1≤i,j≤k) is a matrix of order nd×nd. Here,(4.4)Hi,j={L̃2+Ṽ̃1,i=j=1,M̃2,1≤i≤k-1,j=i+1,L̃2,1<i=j<k,Ṽ̃i,1<i<k,j=1,R1L̃1+Ṽ̃k,i=k,j=1,Ri-1M̃1+RiL̃1,i=k,1<j<k,Rk-1M̃1+RkL̃1+L̃2,i=j=k.We can recursively defineUk=Ũk=Hk-1,k(I-Hk,k)-1, Ui=Hi-1,i(I-Hi,i-Ũi+1Hk,i)-1 for 2≤i≤k-1, Ũi=UiŨi+1 for 2≤i≤k-1, Ũ̃2=U2, Ũ̃i=Ũ̃i-1Ui for 3≤i≤k, and U1=H1,1+∑i=2kŨ̃iHi,1.Let us define another invariant probability vectorũ=[ũ1,ũ2,ũ3,…,ũk] of the stochastic matrix H where ũi(1≤i≤k) is a row vector having nd scalar elements. Here, ũ is computed from ũ1U1=ũ1, ũ11nd=1 and ũi=ũi-1Ui for 2≤i≤k. Next a positive real number u is computed as u=ũ(η+ζ) where η=E(Ink-B0)-11nk and ζ=(Inkd-R)-11nkd. η and ζ are both of size nkd×1. The column vectors η and ζ can be expressed as η=[η1,η2,η3,…,ηk-1,ηk] and ζ=[ζ1,ζ2,ζ3,…,ζk-1,ζk], respectively. Both ηi and ζi(1≤i≤k) are of size nd×1. ηk=Ṽ11n and ηj=ηj+1+Ṽk-j+11n for (1≤j≤k-1). There are only (n-1) nonzero elements in ηk and n nonzero elements in ηj for 1≤j≤k-1. On the other hand, ζi=1nd for 1≤i≤k-1 and ζk=(∑j=1k-1Φj+Υ)1nd.x1 can be computed as x1=u-1ũ using O(n3kd3) multiplications instead of O(n3k3d3) multiplications using the special structure of H. Let us partition x0 and x1 as x0=[κ1,κ2,κ3,…,κk] and x1=[ν1,ν2,ν3,…,νk], respectively, for efficient computation of x0. Here, κi and νi for 1≤i≤k are of size 1×n and 1×nd, respectively. κi(1≤i≤k) is computed as κi=∑j=1iνi-j+1Ṽj and these vector-matrix multiplications can be made efficient by considering only nonzero rows of Ṽj(1≤j≤k).xi(i≥2) can be computed as xi=xi-1R for (i≥2) using n(d-1)×nkd=n2kd(d-1) multiplications instead of n2k2d2 multiplications as R has only n(d-1) nonzero rows.


## 5. Waiting-Time Distribution ofGkI/D/1 System

The waiting time distribution of an arbitrary customer inGIk/D/1 system is the same as the waiting time distribution of a customer in GI/D/k system. Let w(r) be the probability that a customer (i.e., the kth customer of a supercustomer group) has to wait r units of time. w(r) is computed after obtaining the probability vectors yi's where yi(i≥0) is the stationary probability vector that a supercustomer sees i supercustomers ahead of him. y0 can be expressed as y0=[y0,1,y0,2,…,y0,nk]. Similarly, yi(i≥1) can be expressed as yi=[yi,1,yi,2,…,yi,nk] and yi,v can be expressed as yi,v=[yi,v,1,yi,v,2,…,yi,v,d].Two approaches for computing the distribution of waiting time are presented in Sections5.1 and 5.2, respectively. The first method uses the approach for phase type service time distribution and this approach is used in the analysis of Alfa [3] and Chaudhry et al. [7]. In the second method we develop a computationally simpler approach. This is a simplified version of Method 1. It simplifies the computations of the probability vectors yi's for i≥0 and w(r), the probability mass function of the distribution of waiting time.


### 5.1. Method 1

The expressions for the stationary distribution of an entering supercustomer to find different number of supercustomers ahead of him in the system can be written as(5.1)yi={1λ*[x0C1+x1(C1⊗s)],i=0,1λ*[xi(C1⊗S)+xi+1(C1⊗(sβ))],i≥1.In this method, the probability that the waiting time isr units of time is computed as(5.2)w(r)={y01nk,r=0,∑i=1ryi(1nk⊗Id)Ω(i)(r)1d,r≥1,
where Ω(i)(r) is a square matrix of order d that represents the probability distribution of r unit of work in system at the arrival of a supercustomer who finds that i supercustomers are ahead of him. For example, the element (Ω(i)(r))u,v represents the probability that the service of the arbitrary supercustomer, who arrives to find i supercustomers in the system with remaining workload of r units begins in phase v(1≤v≤d), given that service of the supercustomer who was in service at the arrival of the arbitrary supercustomer was in phase u(1≤u≤d). The matrix Ω(i)(r) is computed in the following way:
(5.3)Ω(i)(r)={Id,i=r=0,(sβ)r,1≤i=r,0d,d,i=0,r≥1,0d,d,i≥1,r=0,(S)r-1(sβ),i=1,r≥1,(sβ)Ω(i-1)(r-1)+SΩ(i)(r-1),2≤i≤r.


### 5.2. Method 2

This method simplifies the computation ofy0 and yi(i≥1). Moreover, this method does not need to compute the matrix Ω(i)(r)(i,r≥0) which simplifies the computation of the probability mass function of the distribution of waiting time.The special structures ofC1, s, S, and β can be exploited to reduce the number of operations to compute y0 and yi(i≥1). Only the first columns of the matrices C1 and C1⊗s are nonzero. Therefore, the first element of y0 is nonzero and the remaining (nk-1) elements of y0 are equal to zero. This phenomenon is quite intuitive as the arriving supercustomer does not find anybody ahead of him, he can start his service in phase 1 instantly and the arrival process for next supercustomer starts at phase 1 of the nk phases. The first element y0,1 of y0 is computed as(5.4)y0,1=1λ*∑j=1n(x0,(k-1)n+jtj+x1,(k-1)n+j,dtj).(d-1) columns of (C1⊗S) are nonzero which are column 2 to column d. On the other hand, only the first column of (C1⊗(sβ)) is nonzero. Therefore, the first d elements of yi(i≥1) are nonzero. It is quite intuitive that the first d elements of yi(i≥1) are nonzero (i.e., yi,1 is nonzero and yi,v=01,d(2≤v≤nk)) as the arriving supercustomer finds i supercustomers ahead of him, the arrival process of next supercustomer starts at phase 1 of nk phases and the ongoing service is in the start of phase j(1≤j≤d). The d nonzero elements yi,1,j's (i≥1,1≤j≤d) of the vector yi can be computed using(5.5)yi,1,j={1λ*∑l=1nxi+1,(k-1)n+l,dtl,j=1,1λ*∑l=1nxi,(k-1)n+l,j-1tl,2≤j≤d.The special structures of1nk⊗Id and yi for i≥1 result in(5.6)yi(1nk⊗Id)=yi,1.If a supercustomer seesi supercustomers ahead of him then he has to wait or id+1,id+2,…,id+d amount of time for d>1 (or i units of time for d=1) before getting service.First we consider the case ofd=1. In this case, S=[0], s=[1], and β=[1]. For d=1, (5.3) can be written as(5.7)Ω(i)(r)={1,i=r=0,1,1≤i=r,0,i=0,r≥1,0,i=1,r≥1,Ω(i-1)(r-1),2≤i≤r.Fori<r we can assume i+i1=r and i1>0. In this case, Ω(i)(r)=Ω(i-1)(r-1)=⋯=Ω(0)(i1)=0.Fori=r, we have Ω(i)(r)=1 from (5.7).Fori>r we can assume i=r+i1 and i1>0. In this case, Ω(i)(r)=Ω(i-1)(r-1)=⋯=Ω(i1)(0)=0. Therefore, we can express for d=1,(5.8)Ω(i)(r)1d={1,i=r,0,otherwise.Now, we consider the cased≥2. Let us define ej(1≤j≤d) as the jth column of an identity matrix of order d. It is found that for i≥1 and d≥2(5.9)Ω(i)((i-1)d-j)=0d,dforj=0,1,2,…,(i-1)d,Ω(i)((i-1)d+j)=Sj-1(sβ)=[ed-j+10d,d-1]forj=1,2,…,d,Ω(i)(id+j)=0d,dforj≥1,
which are proved by induction in [30]. If the matrix Ω(i)(r) is nonzero then it contains nonzero elements in the first column. The nonzero first column in nonzero Ω(i)(r) is due to elapsed time representation of service time where service starts at phase 1.Now,[ed-j+10d,d-1]1d=[ed-j+1]. Therefore, we can express for d≥2,(5.10)Ω(i)(r)1d={[ed-(rmodd)],i=⌈rd⌉,0d,1,otherwise.Alfa [2, 3] pointed out that (Ω(i)(r))u,v=0 for i≥1,2≤v≤d. However, he did not provide any explicit expression for Ω(i)(r)1d and his method needs to compute the first column of the matrix Ω(i)(r)(i≥1) recursively as our method does not need to compute. Use of (5.10) to compute Ω(i)(r)1d can significantly reduce the complexity of computing the distribution of waiting time for any multiserver queueing system in discrete time with deterministic service time such as BMAP/D/k system of Alfa [2] and MAP/D/k systems of Alfa [3] and Chaudhry et al. [7]. Equation (5.10) reduces the usage of memory as there is no need to compute and store the matrix Ω(i)(r) for different values of i and r.Now, if an arriving customer has to waitr units of time (r≥0) then there are ⌈r/d⌉ supercustomers ahead of him and the customer, who is in service in the server, is in phase (d-(rmodd)) of service time. Therefore, the probability of a customer waiting r (r≥0) units of time can be expressed by using value of Ω(i)(r)1d of (5.8) and (5.10) in  (5.2)(5.11)w(r)={y0,1,r=0,yr,1,1,r≥1,d=1,y⌈r/d⌉,1,d-(rmodd),r≥1,d≥2.


## 5.1. Method 1

The expressions for the stationary distribution of an entering supercustomer to find different number of supercustomers ahead of him in the system can be written as(5.1)yi={1λ*[x0C1+x1(C1⊗s)],i=0,1λ*[xi(C1⊗S)+xi+1(C1⊗(sβ))],i≥1.In this method, the probability that the waiting time isr units of time is computed as(5.2)w(r)={y01nk,r=0,∑i=1ryi(1nk⊗Id)Ω(i)(r)1d,r≥1,
where Ω(i)(r) is a square matrix of order d that represents the probability distribution of r unit of work in system at the arrival of a supercustomer who finds that i supercustomers are ahead of him. For example, the element (Ω(i)(r))u,v represents the probability that the service of the arbitrary supercustomer, who arrives to find i supercustomers in the system with remaining workload of r units begins in phase v(1≤v≤d), given that service of the supercustomer who was in service at the arrival of the arbitrary supercustomer was in phase u(1≤u≤d). The matrix Ω(i)(r) is computed in the following way:
(5.3)Ω(i)(r)={Id,i=r=0,(sβ)r,1≤i=r,0d,d,i=0,r≥1,0d,d,i≥1,r=0,(S)r-1(sβ),i=1,r≥1,(sβ)Ω(i-1)(r-1)+SΩ(i)(r-1),2≤i≤r.


## 5.2. Method 2

This method simplifies the computation ofy0 and yi(i≥1). Moreover, this method does not need to compute the matrix Ω(i)(r)(i,r≥0) which simplifies the computation of the probability mass function of the distribution of waiting time.The special structures ofC1, s, S, and β can be exploited to reduce the number of operations to compute y0 and yi(i≥1). Only the first columns of the matrices C1 and C1⊗s are nonzero. Therefore, the first element of y0 is nonzero and the remaining (nk-1) elements of y0 are equal to zero. This phenomenon is quite intuitive as the arriving supercustomer does not find anybody ahead of him, he can start his service in phase 1 instantly and the arrival process for next supercustomer starts at phase 1 of the nk phases. The first element y0,1 of y0 is computed as(5.4)y0,1=1λ*∑j=1n(x0,(k-1)n+jtj+x1,(k-1)n+j,dtj).(d-1) columns of (C1⊗S) are nonzero which are column 2 to column d. On the other hand, only the first column of (C1⊗(sβ)) is nonzero. Therefore, the first d elements of yi(i≥1) are nonzero. It is quite intuitive that the first d elements of yi(i≥1) are nonzero (i.e., yi,1 is nonzero and yi,v=01,d(2≤v≤nk)) as the arriving supercustomer finds i supercustomers ahead of him, the arrival process of next supercustomer starts at phase 1 of nk phases and the ongoing service is in the start of phase j(1≤j≤d). The d nonzero elements yi,1,j's (i≥1,1≤j≤d) of the vector yi can be computed using(5.5)yi,1,j={1λ*∑l=1nxi+1,(k-1)n+l,dtl,j=1,1λ*∑l=1nxi,(k-1)n+l,j-1tl,2≤j≤d.The special structures of1nk⊗Id and yi for i≥1 result in(5.6)yi(1nk⊗Id)=yi,1.If a supercustomer seesi supercustomers ahead of him then he has to wait or id+1,id+2,…,id+d amount of time for d>1 (or i units of time for d=1) before getting service.First we consider the case ofd=1. In this case, S=[0], s=[1], and β=[1]. For d=1, (5.3) can be written as(5.7)Ω(i)(r)={1,i=r=0,1,1≤i=r,0,i=0,r≥1,0,i=1,r≥1,Ω(i-1)(r-1),2≤i≤r.Fori<r we can assume i+i1=r and i1>0. In this case, Ω(i)(r)=Ω(i-1)(r-1)=⋯=Ω(0)(i1)=0.Fori=r, we have Ω(i)(r)=1 from (5.7).Fori>r we can assume i=r+i1 and i1>0. In this case, Ω(i)(r)=Ω(i-1)(r-1)=⋯=Ω(i1)(0)=0. Therefore, we can express for d=1,(5.8)Ω(i)(r)1d={1,i=r,0,otherwise.Now, we consider the cased≥2. Let us define ej(1≤j≤d) as the jth column of an identity matrix of order d. It is found that for i≥1 and d≥2(5.9)Ω(i)((i-1)d-j)=0d,dforj=0,1,2,…,(i-1)d,Ω(i)((i-1)d+j)=Sj-1(sβ)=[ed-j+10d,d-1]forj=1,2,…,d,Ω(i)(id+j)=0d,dforj≥1,
which are proved by induction in [30]. If the matrix Ω(i)(r) is nonzero then it contains nonzero elements in the first column. The nonzero first column in nonzero Ω(i)(r) is due to elapsed time representation of service time where service starts at phase 1.Now,[ed-j+10d,d-1]1d=[ed-j+1]. Therefore, we can express for d≥2,(5.10)Ω(i)(r)1d={[ed-(rmodd)],i=⌈rd⌉,0d,1,otherwise.Alfa [2, 3] pointed out that (Ω(i)(r))u,v=0 for i≥1,2≤v≤d. However, he did not provide any explicit expression for Ω(i)(r)1d and his method needs to compute the first column of the matrix Ω(i)(r)(i≥1) recursively as our method does not need to compute. Use of (5.10) to compute Ω(i)(r)1d can significantly reduce the complexity of computing the distribution of waiting time for any multiserver queueing system in discrete time with deterministic service time such as BMAP/D/k system of Alfa [2] and MAP/D/k systems of Alfa [3] and Chaudhry et al. [7]. Equation (5.10) reduces the usage of memory as there is no need to compute and store the matrix Ω(i)(r) for different values of i and r.Now, if an arriving customer has to waitr units of time (r≥0) then there are ⌈r/d⌉ supercustomers ahead of him and the customer, who is in service in the server, is in phase (d-(rmodd)) of service time. Therefore, the probability of a customer waiting r (r≥0) units of time can be expressed by using value of Ω(i)(r)1d of (5.8) and (5.10) in  (5.2)(5.11)w(r)={y0,1,r=0,yr,1,1,r≥1,d=1,y⌈r/d⌉,1,d-(rmodd),r≥1,d≥2.


## 6. Numerical Examples

This section compares our proposed method of computing the matrixR with natural, traditional, and U-based methods and Logarithmic Reduction technique [26] for the matrices R and G. Three different matrices R, G, and U are computed in each of these methods. The values of the convergence constant ε used are 10-6, 10-7, 10-8, and 10-9. Three identical servers (i.e., k=3) are considered with duration of service as eight time units (i.e., d=8) and the support for interarrival times is assumed to be four (i.e., n=4). Three different probability vectors [0.1,0.2,0.2,0.5], [0.1,0.2,0.3,0.4], and [0.1,0.3,0.4,0.2] are used for interarrival times for numerical examples. A program is written in C programming language for different methods to compute the matrices R, G, and U. The program was executed in SunOS 5.10 [31] on an i386 processor which operates at 2660 MHz and has an i387 compatible floating point processor. The number of iterations and time in seconds is recorded in Tables 1–6 for each method from the execution of the program. “LR” is used to denote logarithmic reduction technique in those tables.Table 1
Number of iterations required for different methods for interarrival times having probability vectora¯=[0.10.20.20.5].
NaturalTraditionalU-basedNaturalTraditionalU-basedLRLRStructureεRRRGGGRGR10-61573551944053215710-71864162224653218610-82144762505253221410-924353627958632243Table 2
Time in seconds required for different methods to execute for interarrival times having probability vectora¯=[0.10.20.20.5].
NaturalTraditionalU-basedNaturalTraditionalU-basedLRLRStructureεRRRGGGRGR10-63.18251.11150.20813.87661.26420.20750.34970.25790.072410-73.80691.30130.24514.49461.41590.20700.35300.27510.077610-84.34521.46410.24545.04301.60420.20700.36540.25900.084410-94.95891.66190.26335.69281.78400.24600.35100.25910.0895Table 3
Number of iterations required for different methods for interarrival times having probability vectora¯=[0.1,0.2,0.3,0.4].
NaturalTraditionalU-basedNaturalTraditionalU-basedLRLRStructureεRRRGGGRGR10-62074862535453220710-72455662916253224510-82836473287063228310-932272736679632322Table 4
Time in seconds required for different methods to execute for interarrival times having probability vectora¯=[0.1,0.2,0.3,0.4].
NaturalTraditionalU-basedNaturalTraditionalU-basedLRLRStructureεRRRGGGRGR10-64.24061.48660.25005.11001.65550.20650.35340.25910.083110-75.02941.72830.24556.15091.93750.21990.35340.25860.090410-85.75461.95520.28096.58592.11620.24280.35480.25640.100310-96.50222.19190.27917.31622.37850.24780.35330.25860.1077Table 5
Number of iterations required for different methods for interarrival times having probability vectora¯=[0.1,0.3,0.4,0.2].
NaturalTraditionalU-basedNaturalTraditionalU-basedLRLRStructureεRRRGGGRGR10-618645452522726142254186410-722866562826617212654228510-827077663230518282964270610-9312887636344193433653127Table 6
Time in seconds required for different methods to execute for interarrival times having probability vectora¯=[0.1,0.3,0.4,0.2].
NaturalTraditionalU-basedNaturalTraditionalU-basedLRLRStructureεRRRGGGRGR10-637.558316.24500.911645.214418.23920.79280.53950.44610.457910-746.307619.61711.027853.372221.47090.94530.54570.45850.552610-854.596422.90921.156261.047424.56551.03680.63340.44790.648810-963.937026.91601.303570.162728.38901.20200.63460.54530.7593Tables1 and 2 show the number of iterations and time required to compute the matrices R, G, and U for the interarrival times having probability vector [0.1,0.2,0.2,0.5]. The means of interarrival time and arrival rate for this probability vector are 3.1 time slots and 0.3226/time slot, respectively, and ρ=0.8602. Tables 3 and 4 show the number of iterations and time required for the interarrival times having probability vector [0.1,0.2,0.3,0.4]. The means of interarrival time and arrival rate for the probability vector [0.1,0.2,0.3,0.4] are 3.0 time slots and 0.3333/time slot, respectively, and ρ=0.8889. Our methods of computing the matrix R are found to require least time for these two examples among the methods used here though logarithmic reduction technique and U-based methods require less number of iterations. Tables 5 and 6 show the number of iterations and time required for the interarrival times having probability vector [0.1,0.3,0.4,0.2]. The means of interarrival time and arrival rate for this probability vector [0.1,0.3,0.4,0.2] of interarrival times are 2.7 time slots and 0.3704/time slot, respectively, and ρ=0.9877. Table 6 shows that for the interarrival times probability vector [0.1,0.3,0.4,0.2], Logarithmic Reduction techniques for the matrices G and R require less time to converge than our method of exploiting the structure of the matrices A2, A1, and A0. The relative advantage of our method of exploiting the structures to Logarithmic Reduction techniques decreases with respect to time with increasing value of ρ.


## 7. Conclusion

In this paper, a class of GI/D/k systems in discrete time is analyzed by converting it to a single server queue problem with a convoluted arrival process. Then the stationary distribution of the length of the queue of the transformed system is analyzed using QBD approach. The special structures of the matrices make the computation of the matrix R and the distribution of queue length of the transformed system efficient. Numerical examples show that our proposed method for computing the matrix R can be as efficient as quadratic convergent algorithms such as logarithmic reduction technique. It is also shown that for deterministic service time the waiting time distribution does not need any cumbersome computation like phase type service time distribution.



---

*Source: 716364-2009-11-05.xml*