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

Abstract

A class of discrete time GI/D/𝑘 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.

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/𝑘 system with finite support on interarrival times.

To the best of our knowledge GI/D/𝑘 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/𝑘 systems in discrete time with focus on computational efficiency. The class of GI/D/𝑘 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/𝑐 system by converting the system to an E_𝑐/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/𝑟𝑘 queue with FIFO into E_𝑘/D/𝑟 queues with FIFO for analyzing waiting time distribution. Franx [8] developed expression for the waiting time distribution of the M/D/𝑐 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/𝑁 system in continuous time using spectral technique. Takine [18] analyzed a MAP/D/𝑠 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/𝑐 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/𝑘 system in discrete time. Later Alfa [3] gave a simple and more efficient computational scheme for the MAP/D/𝑘 system in discrete time. Alfa [2] also carried out algorithmic analysis for a discrete time BMAP/D/𝑘 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/𝑐 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], GI^𝑋/D/𝑐 queue by Chaudhry and Kim [5], M^𝑋/D/𝑐 queue by Chaudhry et al. [6] and Franx [9], Ph/D/𝑐 queues by van Hoorn [16], and an E_𝑘/D/𝑟 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/𝑘 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/𝑘 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 𝐑, (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 𝐑 is reduced by decreasing the number of multiplications in each iteration from 𝑂(𝑛³𝑘³𝑑³) to 𝑂(𝑛³𝑘𝑑²) while requiring the same number of iteration as the natural iteration. Here, 𝑛 is the support of general interarrival times, 𝑘 is the number of servers, and 𝑑 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. Section 2 introduces the GI/D/𝑘 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 𝐑. 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 𝐑. Section 7 concludes the paper.

2. The GI/D/𝑘 System

The class of GI/D/𝑘 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, 𝐴, of this multiserver system are general, independent, and identically distributed (iid) with distribution 𝑎_𝑣=Pr{𝐴=𝑣}, 𝑣=1,2,….

We let 𝜆=1/𝐸[𝐴] 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 GI^𝑋/G/1 systems by Alfa [24]. General independent interarrival times can be represented as an absorbing Markov chain with a transition matrix 𝐏_𝐀=[_𝟎^𝐓𝐭_1,𝑛1] where 𝐓 is a square matrix of order 𝑛, 𝐭 is a column vector of order 𝑛, 𝐭=𝟏_𝑛−𝐓𝟏_𝑛, 𝟎_𝑖,𝑗 is an 𝑖×𝑗 matrix of zeros, and 𝟏_𝑖 is a column vector of ones of order 𝑖. The case of infinite interarrival times can be captured for 𝑛=∞. However, the interarrival times are finite if data is collected from practical systems. Therefore, general independent interarrival times with finite support (𝑛<∞) are considered in this paper. The probabilities of interarrival times can be represented by a vector 𝐚=[𝑎₁,𝑎₂,…,𝑎_𝑛]. Here, 𝐚𝟏_𝑛=1 and 𝑎₀=0. This general arrival process can be defined by phase type distribution (𝛼,𝐓) of dimension 𝑛 in elapsed time representation. The parameters 𝛼, 𝐓, and 𝐭 are given as 𝛼=[1,0,0…,0], 𝐓_𝑖,𝑗=(1−∑^𝑖_𝑗=0𝑎_𝑗)/(1−∑^𝑖−1_𝑣=0𝑎_𝑣)=𝑎_𝑖, for 1≤𝑖<𝑛, 𝑗=𝑖+1, 𝐓_𝑖,𝑗=0, for all 𝑗≠𝑖+1 and 𝐭=[𝑎₁𝑎₂⋯𝑎_𝑛−11]^′ where 𝑎_𝑖=1−𝑎_𝑖 for 1≤𝑖≤𝑛−1 and 𝐖^′ is the transpose of matrix 𝐖.

The general independent arrival process can be represented by two square matrices 𝐃₀ and 𝐃₁ of order 𝑛 where (𝐃_𝑙)_𝑖𝑗(𝑙=0,1and1≤𝑖,𝑗≤𝑛) represents a transition from phase 𝑖 to phase 𝑗 with 𝑙 arrivals. Now, 𝐃₀ and 𝐃₁ can be represented in terms of phase type distribution (𝛼,𝐓). Here, 𝐃₀=𝐓 and 𝐃₁=𝐭𝛼=[𝐭𝟎_{𝑛,𝑛−1}]. The first column of 𝐃₁ is nonzero and this column is the vector 𝐭. The matrices 𝐃₀ and 𝐃₁ are similar to the matrices of zero or one arrival for MAP distributed interarrival times. However, they are special cases because the matrices 𝐃₀ and 𝐃₁ only capture general and independently distributed interarrival times.

A discrete time Markov chain {(𝑁(𝑡),𝐽(𝑡));𝑡=0,1,2,…;𝑁(𝑡)∈{0,1,2,…};𝐽(𝑡)∈{1,2,3,…,𝑛}} can be considered to represent this general arrival process where 𝑁(𝑡) is the number of arrivals up to and including time 𝑡 and 𝐽(𝑡) is the phase of the next arrival at time 𝑡.

The matrix 𝐃=𝐃₀+𝐃₁ is stochastic. If 𝜋=[𝜋₁,𝜋₂,…,𝜋_𝑛] is the invariant probability vector for arrival in 𝑛 different phases then 𝜋=𝜋𝐃 and 𝜋𝟏_𝑛=1. The arrival rate is given by 𝜆=𝜋𝐃₁𝟏_𝑛=𝜋𝐭.

There are 𝑘 identical servers in a GI/D/𝑘 system. The service times are of constant duration 𝑑 units (𝑑≥1). Service time can be represented as phase type distribution in elapsed time form with representation (𝛽,𝐒) of dimension 𝑑 where 𝛽=[1,0,0,…,0] and 𝐒=[^{𝟎𝑑−1,1𝐈𝑑−1}_{0𝟎1,𝑑−1}], 𝐈_𝑗 is an identity matrix of order 𝑗. Another column vector 𝐬=𝟏_𝑑−𝐒𝟏_𝑑 can be defined for this phase type distribution (𝛽,𝐒).

The condition that ensures the stability of GI/D/𝑘 system is 𝜌<1 where 𝜌=𝜆𝑑/𝑘. 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/𝑘 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 𝑗th (1≤𝑗≤𝑘) server. In this approach if the first customer is assigned to the first server then the 𝑗th server will be assigned with 𝑗, 𝑗+𝑘, 𝑗+2𝑘,…th customers. In this case we only need to study one server, an arbitrary server 𝑗  (1≤𝑗≤𝑘). This arbitrary server sees the arrival of customers to be according to another independent general process which we call _𝑘GI. This _𝑘GI arrival process is a 𝑘-fold convolution of original arrival process to the multiserver system. _𝑘GI can be defined by two square matrices 𝐂₀ and 𝐂₁ of order 𝑛𝑘. The element (𝐂_𝑣)_𝑖𝑗 (𝑣=0,1 and 1≤𝑖,𝑗≤𝑛𝑘) represents transition from phase 𝑖 to phase 𝑗 with 𝑣 arrivals. Both 𝐂₀ and 𝐂₁ can be considered to comprise of 𝑘² square block matrices of order 𝑛. In this representation (𝐂₀)_𝑖,𝑗=𝐃₀ for 𝑖=𝑗,1≤𝑖≤𝑘, (𝐂₀)_𝑖,𝑗=𝐃₁ for 𝑗=𝑖+1,1≤𝑖<𝑘 and (𝐂₀)_𝑖,𝑗=𝟎_𝑛,𝑛 for 1≤𝑖≤𝑘 and 𝑖≠𝑗 or 𝑗≠𝑖+1. Similarly (𝐂₁)_𝑖,𝑗=𝐃₁ for 𝑖=𝑘,𝑗=1 and (𝐂₁)_𝑖,𝑗=𝟎_𝑛,𝑛 in all other cases for 1≤𝑖,𝑗≤𝑘. The first column and the last row of 𝐂₀ are zero as the first column and last row of 𝐃₀ are zero. Only the last 𝑛 rows of 𝐂₁ are nonzero and only the first element of each of these 𝑛 rows is nonzero.

The matrices 𝐂₀ and 𝐂₁ exhibit behavior similar to the matrices 𝐃₀ and 𝐃₁ except that 𝐂_𝑣 (𝑣=0,1) is a new matrix representing the arrival of 𝑣 supercustomers where a supercustomer is the last customer to form a group of 𝑘 customers or the arrival of 𝑣 regular customers to an arbitrary server. On the other hand, 𝐷_𝑣  (𝑣=0,1) represents arrival of 𝑣 regular customers to the GI/D/𝑘 system. Now, this new arrival process _𝑘GI to an arbitrary server can be represented by a discrete time Markov chain {(𝑁(𝑡),𝐽(̂𝑡));𝑡=0,1,2,…;𝑁(𝑡)∈{0,1,2,…};𝐽(̂𝑡)∈{1,2,…,𝑛𝑘}} where 𝑁(𝑡) is the number of supercustomers that have arrived by time 𝑡 and 𝐽(̂𝑡) is the phase of the next arrival at time 𝑡.

The matrix 𝐂=𝐂₀+𝐂₁ is stochastic. If 𝛿, a 1×𝑛𝑘 matrix, is the invariant probability vector for arrival in different 𝑛𝑘 phases then 𝛿=𝛿𝐂 and 𝛿𝟏_𝑛𝑘=1. 𝛿 can also be represented as 𝛿=[𝛿₁,𝛿₂,𝛿₃,…,𝛿_𝑘] where 𝛿_𝑖(1≤𝑖≤𝑘) is of size 1×𝑛. The supercustomer arrival rate is given by 𝜆^∗=𝜆/𝑘=𝛿𝐂₁𝟏_𝑛𝑘=𝛿_𝑘𝐭.

2.2. The New Model I/D/1 System

A model can be developed for a single server queue with arrival _𝑘GI and deterministic service time of 𝑑 units following the techniques of Alfa [3]. This single server queue is called _𝑘GI/D/1 system and the resulting model is of QBD process type for which standard methods can be used [22]. Let us assume that 𝐿_𝑡, 𝐽_𝑡, and 𝑆_𝑡 represent the number of supercustomers, the arrival phase of the next supercustomer, and the elapsed time of the supercustomer in service, respectively, at time 𝑡 in this QBD model. The state space can be represented by {𝐿_𝑡,𝐽_𝑡,𝑆_𝑡} for 𝐿_𝑡≥1 and the state space is {0,𝐽_𝑡} for 𝐿_𝑡=0. The total state space is {(0,𝑗)∪(𝑖,𝑗,𝑙)∶𝑖≥1,1≤𝑗≤𝑛𝑘,1≤𝑙≤𝑑}. The transition matrix for this _𝑘GI/D/1 system is

𝐁𝐏=⁰𝐁₁𝐄𝐀₁𝐀₀𝐀₂𝐀₁𝐀₀𝐀₂𝐀₁𝐀₀,⋱⋱⋱ (2.1) where 𝐁₀=𝐂₀, 𝐁₁=𝐂₁⊗𝛽, 𝐄=𝐂₀⊗𝐬, 𝐀₂=𝐂₀⊗(𝐬𝛽), 𝐀₁=𝐂₀⊗𝐒+𝐂₁⊗(𝐬𝛽), and 𝐀₀=𝐂₁⊗𝐒.

Let 𝐱=[𝐱₀,𝐱₁,𝐱₂,…] be the stationary distribution of 𝐏 where 𝐱=𝐱𝐏, 𝐱₀𝟏_𝑛𝑘+∑^∞_𝑖=1𝐱_𝑖𝟏_𝑛𝑘𝑑=1 and 𝐱_𝑗 (𝑗≥0) is the probability vector representing 𝑗 supercustomers in the _𝑘GI/D/1 system. During the computation of the stationary distribution of queue length 𝐱, a matrix 𝐑 is computed first which is followed by the computations of 𝐱₀ and 𝐱₁. Then 𝐱_𝑖(𝑖≥2) is computed using the matrix-geometric result 𝐱_𝑖=𝐱_𝑖−1𝐑=𝐱₁𝐑^𝑖−1. Here, the matrix 𝐑 is of order 𝑛𝑘𝑑×𝑛𝑘𝑑 which is the minimal nonnegative solution of 𝐑=𝐀₀+𝐑𝐀₁+𝐑²𝐀₂. The entry 𝐑_𝑢,𝑣(1≤𝑢,𝑣≤𝑛𝑘𝑑) of 𝐑 represents the expected number of visits into (𝑚+1,𝑣), starting from (𝑚,𝑢), before the first return to level 𝑚 (𝑚≥1). Here, 𝑚 is the number of supercustomers in the system and 𝑢 and 𝑣 are any of the 𝑛𝑘𝑑 phases of the system.

3. Structures of the Matrices

This section explains the structures of the block-matrices 𝐁₀, 𝐁₁, 𝐄, 𝐀₂, 𝐀₁, and 𝐀₀ and how their structures are exploited for efficient computation of the matrix 𝐑. The structures of these matrices are also fully exploited in computing the stationary distribution of the queue length in Section 4.

Here, 𝐁₀ and 𝐁₁ are matrices of order 𝑛𝑘×𝑛𝑘 and 𝑛𝑘×𝑛𝑘𝑑, respectively. The first column and the last row of the matrix 𝐁₀ are zero. On the other hand, only the last 𝑛 rows of the matrix 𝐁₁ are nonzero and among these last 𝑛 rows only the first column is nonzero which is 𝐭. The matrix 𝐁₁ can be represented as 𝐁₁=[^{𝟎𝑛(𝑘−1),𝑛𝑑𝟎𝑛(𝑘−1),𝑛𝑑(𝑘−1)}𝐋_{𝟎𝑛,𝑛𝑑(𝑘−1)}] where 𝐋=[𝐭𝟎_{𝑛,𝑛𝑑−1}].

The matrix 𝐄 is of order 𝑛𝑘𝑑×𝑛𝑘 which can be represented as 𝐋𝐌𝐋𝐌𝐋𝐌𝐋𝐄=⋱⋱, (3.1) where 𝑎𝐋=¹𝐬𝑎₂𝐬⋱𝑎_𝑛−1𝐬,𝑎𝐌=¹𝐬𝑎₂𝐬⋮𝑎_𝑛−1𝐬𝐬. (3.2)

𝐋=𝐃₀⊗𝐬 is a matrix of order 𝑛𝑑×𝑛. Moreover, 𝐋 has (𝑛−1) nonzero columns and each column contains only one nonzero element. 𝐋 has (𝑛−1) nonzero elements in (𝑗𝑑,𝑗+1)th (1≤𝑗≤𝑛−1) positions. 𝐌=𝐃₁⊗𝐬 is also a matrix of order 𝑛𝑑×𝑛 matrix. Furthermore, only the first column of 𝐌 is nonzero which contains 𝑛 nonzero elements. 𝐌 has nonzero elements in (𝑗𝑑,1)th (1≤𝑗≤𝑛) positions.

The matrices 𝐀₂, 𝐀₁, and 𝐀₀ are of size 𝑛𝑘𝑑×𝑛𝑘𝑑. The matrix 𝐀₂ can be represented as

𝐀₂=𝐋₁𝐌₁𝐋₁𝐌₁𝐋⋱⋱₁𝐌₁𝐋₁, (3.3) where 𝐋₁ and 𝐌₁ are both matrices of order 𝑛𝑑×𝑛𝑑, 𝐋₁=𝐃₀⊗(𝐬𝛽) and 𝐌₁=𝐃₁⊗(𝐬𝛽). Hence,

𝐋₁=𝑎₁𝑎(𝐬𝛽)₂⋱𝑎(𝐬𝛽)_𝑛−1,𝐌(𝐬𝛽)¹=𝑎₁𝑎(𝐬𝛽)₂⋮𝑎(𝐬𝛽)_𝑛−1.(𝐬𝛽)(𝐬𝛽) (3.4)

There are (𝑛−1) nonzero columns in 𝐋₁ and each column contains only one nonzero element. 𝐋₁ has these nonzero elements in (𝑗𝑑,𝑗𝑑+1)th (1≤𝑗≤𝑛−1) positions. Only the first column of 𝐌₁ is nonzero with 𝑛 nonzero elements in (𝑗𝑑,1)th (1≤𝑗≤𝑛) positions.

The matrix 𝐀₁ can be represented as 𝐀₁=𝐋₂𝐌₂𝐋₂𝐌₂𝐋⋱⋱₂𝐌₂𝐌₁𝐋₂, (3.5) where 𝐋₂=𝑎₁𝐒𝑎₂𝐒⋱𝑎_𝑛−1𝐒,𝐌₂=𝑎₁𝐒𝑎₂𝐒⋮𝑎_𝑛−1𝐒𝐒. (3.6)

Here, 𝐋₂ and 𝐌₂ are both matrices of order 𝑛𝑑×𝑛𝑑. 𝐋₂=𝐃₀⊗𝐒 and 𝐌₂=𝐃₁⊗𝐒. There are (𝑛−1)(𝑑−1) nonzero columns, (𝑛−1)(𝑑−1) nonzero rows, and (𝑛−1)(𝑑−1) nonzero elements in 𝐋₂. Moreover, each nonzero row or each nonzero column of 𝐋₂ contains only one nonzero element. On the other hand, there are (𝑑−1) nonzero columns in 𝐌₂ from column 2 to 𝑑 each having 𝑛 nonzero elements.

The matrix 𝐀₀ can be represented as 𝐀₀=[𝐌₂]. 𝐀₀ has only 𝑛(𝑑−1) nonzero rows which are in its last 𝑛𝑑 rows and (𝑑−1) nonzero columns from column 2 to column 𝑑.

3.1. Computation of the Matrix 𝐑

The matrix 𝐑 can be computed using any of the three linearly convergent formulae from Neuts [22]—𝐑(𝑖+1)∶=𝐀₀+𝐑(𝑖)𝐀₁+𝐑²(𝑖)𝐀₂, 𝐑(𝑖+1)∶=(𝐀₀+𝐑²(𝑖)𝐀₂)(𝐈_𝑛𝑘𝑑−𝐀₁)⁻¹, and 𝐑(𝑖+1)∶=𝐀₀(𝐈_𝑛𝑘𝑑−𝐀₁−𝐑(𝑖)𝐀₂)⁻¹ which are known as natural, traditional, and 𝐔-based iterations, respectively. Here, 𝐑(𝑖) is the 𝑖th iteration value of 𝐑, 𝐑(0)∶=𝟎_{𝑛𝑘𝑑,𝑛𝑘𝑑}, 𝐔 is a square matrix of the same order as the matrix 𝐑 and 𝐔=𝐀₁+𝐑𝐀₂. The matrix 𝐔 is the minimal nonnegative solution of 𝐔=𝐀₁+𝐀₀(𝐈−𝐔)⁻¹𝐀₂. The entry 𝐔_𝑢,𝑣(1≤𝑢,𝑣≤𝑛𝑘𝑑) represents the taboo probability of starting from (𝑚,𝑢) for 𝑚≥1 the eventual visit of the Markov chain to level 𝑚 by visiting (𝑚,𝑣) under the taboo of level (𝑚−1). Any of these three methods terminate in the 𝑖th (𝑖≥1) iteration for ‖𝐑(𝑖)−𝐑(𝑖−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, 𝐔-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 𝐔-based iteration is less than the numbers of iterations for other two methods. Alfa and Xue [25] developed efficient inversion technique for 𝐔-based iteration for the matrix 𝐑 in analyzing the queue length of a GI/G/1 system in discrete time.

There are several efficient methods for computing the matrix 𝐑 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 𝐆 which is the minimal nonnegative solution to 𝐆=𝐀₂+𝐀₁𝐆+𝐀₀𝐆². Here, 𝐆 is an 𝑚×𝑚 matrix for 𝑚×𝑚 matrices 𝐀₂, 𝐀₁, and 𝐀₀ where 𝐺_𝑖,𝑗  (1≤𝑖,𝑗≤𝑚) represents the probability that starting from state (𝑁+1,𝑖) the Markov chain eventually visits the level 𝑁  (𝑁≥1) and does so by visiting the state (𝑁,𝑗). The matrices 𝐑, 𝐔, and 𝐆 are related and these relations can be found in [26]. Similar to the computation of the matrix 𝐑, the matrix 𝐆 can be computed using three different iterations—natural, traditional, and 𝐔-based. Generally a zero matrix is used for the initial value for iterations of the matrix 𝐆 (i.e., 𝐆=𝟎_𝑚,𝑚). Latouche and Ramaswami [26] showed that if 𝑆_𝐺(𝜀) iterative steps are needed to converge for a very small positive constant 𝜀 such that ‖𝐆(𝑆_𝐺(𝜀))−𝐆(𝑆_𝐺(𝜀)−1)‖_∞<𝜀 for the (𝑙+1)-st natural iteration 𝐆(𝑙+1)=𝐀₂+𝐀₁𝐆(𝑙)+𝐀₀𝐆²(𝑙), then logarithmic reduction requires not more than log₂𝑆_𝐺(𝜀) iterations and overall complexity of their algorithm is 𝑂(𝑚³log₂𝑆_𝐺(𝜀)). They also developed logarithmic reduction technique for the matrix 𝐑. However, the quadratic convergent algorithms [26–28] are not considered here for computing the matrix 𝐑 as these algorithms involve inversion of matrices which cannot exploit the special structures of the matrices 𝐀₂, 𝐀₁, and 𝐀₀.

The complexity of computing the matrix 𝐑 can be reduced by using the structure of 𝐀₀. The first 𝑛𝑑(𝑘−1) rows are zero in 𝐀₀. Therefore, the first 𝑛𝑑(𝑘−1) rows are zero in 𝐑. Among the last 𝑛𝑑 rows of 𝐀₀, 𝑛 rows are zero which are rows 𝑛𝑑(𝑘−1)+𝑑, 𝑛𝑑(𝑘−1)+2𝑑, 𝑛𝑑(𝑘−1)+3𝑑,…,𝑛𝑘𝑑. The corresponding rows are zero for 𝐑. The matrix 𝐑 can be represented as follows:

𝐑𝐑=₁𝐑₂𝐑₃⋯𝐑_𝑘. (3.7) Here, 𝐑_𝑗(1≤𝑗≤𝑘) is a matrix of order 𝑛𝑑×𝑛𝑑 which contains 𝑑,2𝑑,3𝑑,…,𝑛𝑑th rows as zero rows.

The computations of 𝐑(𝑖)𝐀₂ and 𝐑(𝑖)𝐀₁ require (2𝑘−1) and 2𝑘 block-matrix multiplications, respectively, where each block matrix is of order 𝑛𝑑×𝑛𝑑. Similarly, the computation of 𝐑²(𝑖)𝐀₂ requires multiplication of 𝐑(𝑖) and 𝐑(𝑖)𝐀₂ which involves 𝑘 block-matrix multiplications. Inversions of both the matrices (𝐈_𝑛𝑘𝑑−𝐀₁) and (𝐈_𝑛𝑘𝑑−𝐀₁−𝐑(𝑖)𝐀₂) generate 𝑘² block matrices each of order 𝑛𝑑×𝑛𝑑. Inversion of the matrix (𝐈_𝑛𝑘𝑑−𝐀₁) needs to be computed only once for traditional iteration whereas inversion of the matrix (𝐈_𝑛𝑘𝑑−𝐀₁−𝐑(𝑖)𝐀₂) needs to be computed in each iteration for 𝐔-based iteration. On the other hand, the multiplication of (𝐀₀+𝐑²(𝑖)𝐀₂) and (𝐈_𝑛𝑘𝑑−𝐀₁)⁻¹ in traditional iteration and the multiplication of 𝐀₀ and (𝐈_𝑛𝑘𝑑−𝐀₁−𝐑(𝑖)𝐀₂)⁻¹ in 𝐔-based iteration involve 𝑘² and 𝑘 block-matrix multiplications, respectively, where each block matrix is of order 𝑛𝑑×𝑛𝑑. Each iteration of natural method requires substantially less multiplications than each step of traditional and 𝐔-based methods do. Therefore, we use natural iteration for developing our algorithm for computing the matrix 𝐑 in this paper.

The matrix 𝐑_𝑗(1≤𝑗≤𝑘) of (3.7) can be computed using the natural iteration with the following equation:

𝐑_𝑗𝐌(𝑖+1)∶=²+𝐑₁(𝑖)𝐋₂+𝐑_𝑘(𝑖)𝐌₁+𝐑_𝑘(𝑖)𝐑₁(𝑖)𝐋₁𝐑,𝑗=1,_𝑗−1(𝑖)𝐌₂+𝐑_𝑗(𝑖)𝐋₂+𝐑_𝑘(𝑖)𝐑_𝑗−1(𝑖)𝐌₁+𝐑_𝑘(𝑖)𝐑_𝑗(𝑖)𝐋₁,1<𝑗≤𝑘, (3.8) where 𝐑_𝑗(𝑖)(𝑖≥0,1≤𝑗≤𝑘) is the 𝑖th iteration value of 𝐑_𝑗 with the initial value 𝐑_𝑗(0)∶=𝟎_{𝑛𝑑,𝑛𝑑}(1≤𝑗≤𝑘). The terminating condition for this iteration is same as the terminating condition for natural, traditional, and 𝐔-based methods but needs to consider only 𝑛(𝑑−1) nonzero rows of (𝐑(𝑖)−𝐑(𝑖−1)).

The number of multiplications required to compute 𝐑_𝑗(𝑖)𝐋₂(1≤𝑗≤𝑘) is 𝑛(𝑛−1)(𝑑−1)². Similarly, the number of multiplications required to compute both 𝐑_𝑗(𝑖)𝐋₁ and 𝐑_𝑗(𝑖)𝐌₁ for 1≤𝑗≤𝑘 is 𝑛(𝑛−1)(𝑑−1). In the same way, the number of multiplications required to compute 𝐑_𝑗−1(𝑖)𝐌₂(1<𝑗≤𝑘) is 𝑛(𝑛−1)(𝑑−1)². On the other hand, the numbers of multiplications required to compute 𝐑_𝑘(𝑖)𝐑_𝑗−1(𝑖)𝐌₁(1<𝑗≤𝑘) by block-matrix multiplication of 𝐑_𝑘(𝑖) and 𝐑_𝑗−1(𝑖)𝐌₁(1<𝑗≤𝑘) and 𝐑_𝑘(𝑖)𝐑_𝑗(𝑖)𝐋₁(1≤𝑗≤𝑘) by block-matrix multiplication of 𝐑_𝑘(𝑖) and 𝐑_𝑗(𝑖)𝐋₁(1≤𝑗≤𝑘) are 𝑛²(𝑑−1)² and 𝑛²(𝑛−1)(𝑑−1)², respectively. Hence, the numbers of multiplications required for computing 𝐑₁(𝑖+1) from 𝐑₁(𝑖) and 𝐑_𝑗(𝑖+1)(2≤𝑗≤𝑘) from 𝐑_𝑗(𝑖) in one iteration are 𝑛(𝑛−1)(𝑑−1)(𝑑+1)+𝑛²(𝑛−1)(𝑑−1)² and 𝑛³(𝑑−1)²+2𝑛𝑑(𝑛−1)(𝑑−1), respectively. Therefore, the total number of multiplications required to compute all 𝑘 block matrices 𝐑_𝑗(𝑖+1)'s (1≤𝑗≤𝑘) of the matrix 𝐑(𝑖+1) in one iteration is 𝑛²(𝑛−1)(𝑑−1)²+𝑛(𝑛−1)(𝑑−1)(𝑑+1)+(𝑘−1)(𝑛³(𝑑−1)²+2𝑛𝑑(𝑛−1)(𝑑−1))=𝑂(𝑛³𝑘𝑑²). By exploiting the special structures of the matrices the computational complexity is decreased in each iteration from 𝑂(𝑛³𝑘³𝑑³) to 𝑂(𝑛³𝑘𝑑²) which is substantial reduction in computations. Moreover, the memory requirement for our method is also less than that of natural, traditional, 𝐔-based methods and quadratic convergent algorithms.

4. Stationary Distribution of P for I/D/1 System

Two methods for computing stationary distribution 𝐱 for the Markov chain represented by the matrix 𝐏 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 vector 𝐳=[𝐳₀,𝐳₁] needs to be computed for a stochastic matrix 𝐏=[^𝐁0𝐁1_{𝐄𝐀1+𝐑𝐀2}]. The row vector 𝐳 is computed from 𝐳𝐏=𝐳 and 𝐳₀𝟏_𝑛𝑘+𝐳₁𝟏_𝑛𝑘𝑑=1. Then 𝐱₀ and 𝐱₁ can be expressed as 𝐱₀=𝑣⁻¹𝐳₀ and 𝐱₁=𝑣⁻¹𝐳₁ where 𝑣=𝐳₀𝟏_𝑛𝑘+𝐳₁(𝐈_𝑛𝑘𝑑−𝐑)⁻¹𝟏_𝑛𝑘𝑑.

In the second method, we need to compute 𝐡, the invariant probability vector of a stochastic matrix 𝐇 defined as 𝐇=𝐀₁+𝑅𝐀₂+𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹𝐁₁. Next a positive real number ℎ is computed as follows ℎ=𝐡(𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹𝟏_𝑛𝑘+(𝐈_𝑛𝑘𝑑−𝐑)⁻¹𝟏_𝑛𝑘𝑑). On the other hand, 𝐡 is computed using the relations 𝐡=𝐡𝐇 and 𝐡𝟏_𝑛𝑘𝑑=1. 𝐱₀ and 𝐱₁ can be computed as 𝐱₁=ℎ⁻¹𝐡, 𝐱₀=𝐱₁𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹.

4.1. Exploiting Structures for Stationary Distribution

The stationary distribution of 𝐱 can be computed by first computing the probability vector 𝐱₁ which can be made efficient by exploiting the special structures of the matrices 𝐇, 𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹𝟏_𝑛𝑘 and (𝐈_𝑛𝑘𝑑−𝐑)⁻¹𝟏_𝑛𝑘𝑑. The special structures of the matrices (𝐈_𝑛𝑘𝑑−𝐑)⁻¹, (𝐈_𝑛𝑘𝑑−𝐑)⁻¹𝟏_𝑛𝑘𝑑, (𝐈_𝑛𝑘−𝐁₀)⁻¹, 𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹, 𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹𝟏_𝑛𝑘, 𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹𝐁₁, and 𝐇 are due to the special structures of the matrices 𝐁₀, 𝐁₁, 𝐄, and 𝐑. Here, how the computations of these matrices can be carried out efficiently is explained.

Computing (𝐈_𝑛𝑘𝑑−𝐑)⁻¹=[^{𝐈𝑛𝑑(𝑘−1)𝟎𝑛𝑑(𝑘−1),𝑛𝑑}_Φ𝚼] requires 𝑂(𝑛³𝑘𝑑³) multiplications by using block-matrix inversion [29] where 𝚼=𝐑_𝑘⁻¹, Φ=𝚼Ψ, Ψ=[𝐑₁,𝐑₂,𝐑₃,…,𝐑_𝑘−1], and 𝐑_𝑘=𝐈_𝑛𝑑−𝐑_𝑘. Φ can be expressed as Φ=[Φ₁,Φ₂,Φ₃,…,Φ_𝑘−1] where Φ_𝑗=𝚼𝐑_𝑗(1≤𝑗≤𝑘−1).

(𝐈_𝑛𝑘−𝐁₀)⁻¹ can be represented as

(𝐈_𝑛𝑘−𝐁₀)⁻¹=𝐕₁𝐕₂⋯𝐕_𝑘−1𝐕_𝑘𝐕₁⋯𝐕_𝑘−2𝐕_𝑘−1𝐕⋱⋮⋮₁𝐕₂𝐕₁, (4.1) where 𝐕₁=𝐃₀⁻¹, 𝐕_𝑖=(𝐃₀⁻¹𝐃₁)^𝑖−1𝐃₀⁻¹(1<𝑖≤𝑘), and 𝐃₀=𝐈_𝑛−𝐃₀. The computation of (𝐈_𝑛𝑘−𝐁₀)⁻¹ requires 𝑂(𝑛²𝑘+𝑛𝑘²) instead of 𝑂(𝑛³𝑘³) multiplications by utilizing the special structure of (𝐈_𝑛𝑘−𝐁₀).

The matrices 𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹ and 𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹𝐁₁ can be represented as

𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹=𝐕₁𝐕₂…𝐕_𝑘−1𝐕_𝑘𝐕₁…𝐕_𝑘−2𝐕_𝑘−1𝐕⋱⋮⋮₁𝐕₂𝐕₁,𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹𝐁₁=𝐕₁𝐕₂⋮𝐕_𝑘−1𝐕_𝑘, (4.2) where 𝐕₁=𝐋𝐕₁, 𝐕_𝑖=𝐋𝐕_𝑖+𝐌𝐕_𝑖−1, for 1<𝑖≤𝑘, and 𝐕_𝑖=𝐕_{𝑘−𝑖+1}𝐋(1≤𝑖≤𝑘). Here, 𝐕_𝑖 and 𝐕_𝑖 for 1≤𝑖≤𝑘 are matrices of order 𝑛𝑑×𝑛 and 𝑛𝑑×𝑛𝑑, respectively. Here, multiplying 𝐄 by (𝐈_𝑛𝑘−𝐁₀)⁻¹ and 𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹ by 𝐁₁ both require 𝑂(𝑛²𝑘) multiplications due to utilizing special structures of matrices. If there were no special structures of 𝐄, (𝐈_𝑛𝑘−𝐁₀)⁻¹ and 𝐁₁ then the multiplying 𝐄 by (𝐈_𝑛𝑘−𝐁₀)⁻¹ and 𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹ by 𝐁₁ would require 𝑂(𝑛³𝑘³𝑑) and 𝑂(𝑛³𝑘²𝑑³) multiplications, respectively.

The matrix 𝐇 of can be represented in the following form:

𝐇𝐇=^1,1𝐇_1,2𝐇_2,1𝐇_2,2𝐇_2,3𝐇_3,1𝐇_3,3𝐇_3,4𝐇⋮⋱⋱_𝑘−1,1𝐇_{𝑘−1,𝑘−1}𝐇_{𝑘−1,𝑘}𝐇_𝑘,1𝐇_𝑘,2𝐇_𝑘,3⋯𝐇_{𝑘,𝑘−1}𝐇_𝑘,𝑘, (4.3) where 𝐇_𝑖,𝑗(1≤𝑖,𝑗≤𝑘) is a matrix of order 𝑛𝑑×𝑛𝑑. Here,

𝐇_𝑖,𝑗=𝐋₂+𝐕₁𝐌,𝑖=𝑗=1,₂𝐋,1≤𝑖≤𝑘−1,𝑗=𝑖+1,₂𝐕,1<𝑖=𝑗<𝑘,_𝑖𝐑,1<𝑖<𝑘,𝑗=1,₁𝐋₁+𝐕_𝑘𝐑,𝑖=𝑘,𝑗=1,_𝑖−1𝐌₁+𝐑_𝑖𝐋₁𝐑,𝑖=𝑘,1<𝑗<𝑘,_𝑘−1𝐌₁+𝐑_𝑘𝐋₁+𝐋₂,𝑖=𝑗=𝑘. (4.4)

We can recursively define 𝐔_𝑘=𝐔_𝑘=𝐇_{𝑘−1,𝑘}(𝐈−𝐇_𝑘,𝑘)⁻¹, 𝐔_𝑖=𝐇_{𝑖−1,𝑖}(𝐈−𝐇_𝑖,𝑖−𝐔_𝑖+1𝐇_𝑘,𝑖)⁻¹ for 2≤𝑖≤𝑘−1, 𝐔_𝑖=𝐔_𝑖𝐔_𝑖+1 for 2≤𝑖≤𝑘−1, 𝐔₂=𝐔₂, 𝐔_𝑖=𝐔_𝑖−1𝐔_𝑖 for 3≤𝑖≤𝑘, and 𝐔₁=𝐇_1,1+∑^𝑘_𝑖=2𝐔_𝑖𝐇_𝑖,1.

Let us define another invariant probability vector 𝐮=[𝐮₁,𝐮₂,𝐮₃,…,𝐮_𝑘] of the stochastic matrix 𝐇 where 𝐮_𝑖(1≤𝑖≤𝑘) is a row vector having 𝑛𝑑 scalar elements. Here, 𝐮 is computed from 𝐮₁𝐔₁=𝐮₁, 𝐮₁𝟏_𝑛𝑑=1 and 𝐮_𝑖=𝐮_𝑖−1𝐔_𝑖 for 2≤𝑖≤𝑘. Next a positive real number 𝑢 is computed as 𝑢=𝐮(𝜂+𝜁) where 𝜂=𝐄(𝐈_𝑛𝑘−𝐁₀)⁻¹𝟏_𝑛𝑘 and 𝜁=(𝐈_𝑛𝑘𝑑−𝐑)⁻¹𝟏_𝑛𝑘𝑑. 𝜂 and 𝜁 are both of size 𝑛𝑘𝑑×1. The column vectors 𝜂 and 𝜁 can be expressed as 𝜂=[𝜂₁,𝜂₂,𝜂₃,…,𝜂_𝑘−1,𝜂_𝑘] and 𝜁=[𝜁₁,𝜁₂,𝜁₃,…,𝜁_𝑘−1,𝜁_𝑘], respectively. Both 𝜂_𝑖 and 𝜁_𝑖(1≤𝑖≤𝑘) are of size 𝑛𝑑×1. 𝜂_𝑘=𝐕₁𝟏_𝑛 and 𝜂_𝑗=𝜂_𝑗+1+𝐕_{𝑘−𝑗+1}𝟏_𝑛 for (1≤𝑗≤𝑘−1). There are only (𝑛−1) nonzero elements in 𝜂_𝑘 and 𝑛 nonzero elements in 𝜂_𝑗 for 1≤𝑗≤𝑘−1. On the other hand, 𝜁_𝑖=𝟏_𝑛𝑑 for 1≤𝑖≤𝑘−1 and 𝜁_𝑘=(∑^𝑘−1_𝑗=1Φ_𝑗+𝚼)𝟏_𝑛𝑑.

𝐱₁ can be computed as 𝐱₁=𝑢⁻¹𝐮 using 𝑂(𝑛³𝑘𝑑³) multiplications instead of 𝑂(𝑛³𝑘³𝑑³) multiplications using the special structure of 𝐇. Let us partition 𝐱₀ and 𝐱₁ as 𝐱₀=[𝜅₁,𝜅₂,𝜅₃,…,𝜅_𝑘] and 𝐱₁=[𝜈₁,𝜈₂,𝜈₃,…,𝜈_𝑘], respectively, for efficient computation of 𝐱₀. Here, 𝜅_𝑖 and 𝜈_𝑖 for 1≤𝑖≤𝑘 are of size 1×𝑛 and 1×𝑛𝑑, respectively. 𝜅_𝑖(1≤𝑖≤𝑘) is computed as 𝜅_𝑖=∑^𝑖_𝑗=1𝜈_{𝑖−𝑗+1}𝐕_𝑗 and these vector-matrix multiplications can be made efficient by considering only nonzero rows of 𝐕_𝑗(1≤𝑗≤𝑘).

𝐱_𝑖(𝑖≥2) can be computed as 𝐱_𝑖=𝐱_𝑖−1𝐑 for (𝑖≥2) using 𝑛(𝑑−1)×𝑛𝑘𝑑=𝑛²𝑘𝑑(𝑑−1) multiplications instead of 𝑛²𝑘²𝑑² multiplications as 𝐑 has only 𝑛(𝑑−1) nonzero rows.

5. Waiting-Time Distribution of I/D/1 System

The waiting time distribution of an arbitrary customer in _𝑘GI/D/1 system is the same as the waiting time distribution of a customer in GI/D/𝑘 system. Let 𝑤(𝑟) be the probability that a customer (i.e., the 𝑘th customer of a supercustomer group) has to wait 𝑟 units of time. 𝑤(𝑟) is computed after obtaining the probability vectors 𝐲_𝑖's where 𝐲_𝑖(𝑖≥0) is the stationary probability vector that a supercustomer sees 𝑖 supercustomers ahead of him. 𝐲₀ can be expressed as 𝐲₀=[𝐲_0,1,𝐲_0,2,…,𝐲_0,𝑛𝑘]. Similarly, 𝐲_𝑖(𝑖≥1) can be expressed as 𝐲_𝑖=[𝐲_𝑖,1,𝐲_𝑖,2,…,𝐲_{𝑖,𝑛𝑘}] and 𝐲_𝑖,𝑣 can be expressed as 𝐲_𝑖,𝑣=[𝐲_𝑖,𝑣,1,𝐲_𝑖,𝑣,2,…,𝐲_{𝑖,𝑣,𝑑}].

Two approaches for computing the distribution of waiting time are presented in Sections 5.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 𝐲_𝑖's for 𝑖≥0 and 𝑤(𝑟), 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

𝐲_𝑖=1𝜆^∗[𝐱₀𝐂₁+𝐱₁(𝐂₁1⊗𝐬)],𝑖=0,𝜆^∗[𝐱_𝑖(𝐂₁⊗𝐒)+𝐱_𝑖+1(𝐂₁⊗(𝐬𝛽))],𝑖≥1. (5.1)

In this method, the probability that the waiting time is 𝑟 units of time is computed as

𝐲𝑤(𝑟)=⁰𝟏_𝑛𝑘,𝑟=0,_𝑟𝑖=1𝐲_𝑖(𝟏_𝑛𝑘⊗𝐈_𝑑)Ω^(𝑖)(𝑟)𝟏_𝑑,𝑟≥1, (5.2) where Ω^(𝑖)(𝑟) is a square matrix of order 𝑑 that represents the probability distribution of 𝑟 unit of work in system at the arrival of a supercustomer who finds that 𝑖 supercustomers are ahead of him. For example, the element (Ω^(𝑖)(𝑟))_𝑢,𝑣 represents the probability that the service of the arbitrary supercustomer, who arrives to find 𝑖 supercustomers in the system with remaining workload of 𝑟 units begins in phase 𝑣(1≤𝑣≤𝑑), given that service of the supercustomer who was in service at the arrival of the arbitrary supercustomer was in phase 𝑢(1≤𝑢≤𝑑). The matrix Ω^(𝑖)(𝑟) is computed in the following way: Ω^(𝑖)𝐈(𝑟)=^𝑑,𝑖=𝑟=0,(𝐬𝛽)^𝑟𝟎,1≤𝑖=𝑟,_𝑑,𝑑𝟎,𝑖=0,𝑟≥1,_𝑑,𝑑,𝑖≥1,𝑟=0,(𝐒)^𝑟−1(𝐬𝛽),𝑖=1,𝑟≥1,(𝐬𝛽)Ω^(𝑖−1)(𝑟−1)+𝐒Ω^(𝑖)(𝑟−1),2≤𝑖≤𝑟. (5.3)

5.2. Method 2

This method simplifies the computation of 𝐲₀ and 𝐲_𝑖(𝑖≥1). Moreover, this method does not need to compute the matrix Ω^(𝑖)(𝑟)(𝑖,𝑟≥0) which simplifies the computation of the probability mass function of the distribution of waiting time.

The special structures of 𝐂₁, 𝐬, 𝐒, and 𝛽 can be exploited to reduce the number of operations to compute 𝐲₀ and 𝐲_𝑖(𝑖≥1). Only the first columns of the matrices 𝐂₁ and 𝐂₁⊗𝑠 are nonzero. Therefore, the first element of 𝐲₀ is nonzero and the remaining (𝑛𝑘−1) elements of 𝐲₀ 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 𝑛𝑘 phases. The first element 𝐲_0,1 of 𝐲₀ is computed as

𝐲_0,1=1𝜆^∗𝑛_𝑗=1(𝐱_{0,(𝑘−1)𝑛+𝑗}𝐭_𝑗+𝐱_{1,(𝑘−1)𝑛+𝑗,𝑑}𝐭_𝑗). (5.4)

(𝑑−1) columns of (𝐂₁⊗𝐒) are nonzero which are column 2 to column 𝑑. On the other hand, only the first column of (𝐂₁⊗(𝑠𝛽)) is nonzero. Therefore, the first 𝑑 elements of 𝐲_𝑖(𝑖≥1) are nonzero. It is quite intuitive that the first 𝑑 elements of 𝐲_𝑖(𝑖≥1) are nonzero (i.e., 𝐲_𝑖,1 is nonzero and 𝐲_𝑖,𝑣=𝟎_1,𝑑(2≤𝑣≤𝑛𝑘)) as the arriving supercustomer finds 𝑖 supercustomers ahead of him, the arrival process of next supercustomer starts at phase 1 of 𝑛𝑘 phases and the ongoing service is in the start of phase 𝑗(1≤𝑗≤𝑑). The 𝑑 nonzero elements 𝐲_𝑖,1,𝑗's (𝑖≥1,1≤𝑗≤𝑑) of the vector 𝐲_𝑖 can be computed using

𝐲_𝑖,1,𝑗=1𝜆^∗𝑛_𝑙=1𝐱_{𝑖+1,(𝑘−1)𝑛+𝑙,𝑑}𝐭_𝑙1,𝑗=1,𝜆^∗𝑛_𝑙=1𝐱_{𝑖,(𝑘−1)𝑛+𝑙,𝑗−1}𝐭_𝑙,2≤𝑗≤𝑑. (5.5)

The special structures of 𝟏_𝑛𝑘⊗𝐈_𝑑 and 𝐲_𝑖 for 𝑖≥1 result in

𝐲_𝑖(𝟏_𝑛𝑘⊗𝐈_𝑑)=𝐲_𝑖,1. (5.6)

If a supercustomer sees 𝑖 supercustomers ahead of him then he has to wait or 𝑖𝑑+1,𝑖𝑑+2,…,𝑖𝑑+𝑑 amount of time for 𝑑>1 (or 𝑖 units of time for 𝑑=1) before getting service.

First we consider the case of 𝑑=1. In this case, 𝐒=[0], 𝐬=[1], and 𝛽=[1]. For 𝑑=1, (5.3) can be written as

Ω^(𝑖)Ω(𝑟)=1,𝑖=𝑟=0,1,1≤𝑖=𝑟,0,𝑖=0,𝑟≥1,0,𝑖=1,𝑟≥1,_(𝑖−1)(𝑟−1),2≤𝑖≤𝑟. (5.7)

For 𝑖<𝑟 we can assume 𝑖+𝑖₁=𝑟 and 𝑖₁>0. In this case, Ω^(𝑖)(𝑟)=Ω^(𝑖−1)(𝑟−1)=⋯=Ω⁽⁰⁾(𝑖₁)=0.

For 𝑖=𝑟, we have Ω^(𝑖)(𝑟)=1 from (5.7).

For 𝑖>𝑟 we can assume 𝑖=𝑟+𝑖₁ and 𝑖₁>0. In this case, Ω^(𝑖)(𝑟)=Ω^(𝑖−1)(𝑟−1)=⋯=Ω^(𝑖1)(0)=0. Therefore, we can express for 𝑑=1,

Ω^(𝑖)(𝑟)𝟏_𝑑=1,𝑖=𝑟,0,otherwise. (5.8)

Now, we consider the case 𝑑≥2. Let us define 𝐞_𝑗(1≤𝑗≤𝑑) as the 𝑗th column of an identity matrix of order 𝑑. It is found that for 𝑖≥1 and 𝑑≥2

Ω^(𝑖)((𝑖−1)𝑑−𝑗)=𝟎_𝑑,𝑑forΩ𝑗=0,1,2,…,(𝑖−1)𝑑,_(𝑖)((𝑖−1)𝑑+𝑗)=𝐒^𝑗−1(𝐬𝛽)=[𝐞_{𝑑−𝑗+1}𝟎_{𝑑,𝑑−1}]forΩ𝑗=1,2,…,𝑑,_(𝑖)(𝑖𝑑+𝑗)=𝟎_𝑑,𝑑for𝑗≥1, (5.9) which are proved by induction in [30]. If the matrix Ω^(𝑖)(𝑟) is nonzero then it contains nonzero elements in the first column. The nonzero first column in nonzero Ω^(𝑖)(𝑟) is due to elapsed time representation of service time where service starts at phase 1.

Now, [𝐞_{𝑑−𝑗+1}𝟎_{𝑑,𝑑−1}]𝟏_𝑑=[𝐞_{𝑑−𝑗+1}]. Therefore, we can express for 𝑑≥2,

Ω^(𝑖)(𝑟)𝟏_𝑑=[𝐞_{𝑑−(𝑟mod𝑑)}𝑟],𝑖=⌈𝑑𝟎⌉,_𝑑,1,otherwise. (5.10)

Alfa [2, 3] pointed out that (Ω^(𝑖)(𝑟))_𝑢,𝑣=0 for 𝑖≥1,2≤𝑣≤𝑑. However, he did not provide any explicit expression for Ω^(𝑖)(𝑟)𝟏_𝑑 and his method needs to compute the first column of the matrix Ω^(𝑖)(𝑟)(𝑖≥1) recursively as our method does not need to compute. Use of (5.10) to compute Ω^(𝑖)(𝑟)𝟏_𝑑 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/𝑘 system of Alfa [2] and MAP/D/𝑘 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 Ω^(𝑖)(𝑟) for different values of 𝑖 and 𝑟.

Now, if an arriving customer has to wait 𝑟 units of time (𝑟≥0) then there are ⌈𝑟/𝑑⌉ supercustomers ahead of him and the customer, who is in service in the server, is in phase (𝑑−(𝑟mod𝑑)) of service time. Therefore, the probability of a customer waiting 𝑟 (𝑟≥0) units of time can be expressed by using value of Ω^(𝑖)(𝑟)𝟏_𝑑 of (5.8) and (5.10) in  (5.2)

𝐲𝑤(𝑟)=𝐲^0,1,𝑟=0,_𝑟,1,1𝐲,𝑟≥1,𝑑=1,_{⌈𝑟/𝑑⌉,1,𝑑−(𝑟mod𝑑)},𝑟≥1,𝑑≥2. (5.11)

6. Numerical Examples

This section compares our proposed method of computing the matrix 𝐑 with natural, traditional, and U-based methods and Logarithmic Reduction technique [26] for the matrices 𝐑 and 𝐆. Three different matrices 𝐑, 𝐆, and 𝐔 are computed in each of these methods. The values of the convergence constant 𝜀 used are 10⁻⁶, 10⁻⁷, 10⁻⁸, and 10⁻⁹. Three identical servers (i.e., 𝑘=3) are considered with duration of service as eight time units (i.e., 𝑑=8) and the support for interarrival times is assumed to be four (i.e., 𝑛=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 𝐑, 𝐆, and 𝐔. 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 vector 𝐚=[0.10.20.20.5].

Table 2 

Time in seconds required for different methods to execute for interarrival times having probability vector 𝐚=[0.10.20.20.5].

Table 3 

Number of iterations required for different methods for interarrival times having probability vector 𝐚=[0.1,0.2,0.3,0.4].

Table 4 

Time in seconds required for different methods to execute for interarrival times having probability vector 𝐚=[0.1,0.2,0.3,0.4].

Table 5 

Number of iterations required for different methods for interarrival times having probability vector 𝐚=[0.1,0.3,0.4,0.2].

Table 6 

Time in seconds required for different methods to execute for interarrival times having probability vector 𝐚=[0.1,0.3,0.4,0.2].

Tables 1 and 2 show the number of iterations and time required to compute the matrices 𝐑, 𝐆, and 𝐔 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 𝐑 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 𝐆 and 𝐑 require less time to converge than our method of exploiting the structure of the matrices 𝐀₂, 𝐀₁, and 𝐀₀. 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/𝑘 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 𝐑 and the distribution of queue length of the transformed system efficient. Numerical examples show that our proposed method for computing the matrix 𝐑 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.

Acknowledgment

This research is partially supported by grants from NSERC to A. S. Alfa.