A Discrete Single Server Queue with Markovian Arrivals and Phase Type Group Services

We consider a single-server discrete queucing system in which arrivals occur according to a Markovian arrival process. Service is provided in groups of size no more than M customers. The service times are assumed to follow a discrete phase type distribution, whose representation may depend on the group size. Under a probabilistic service rule, which depends on the number of customers waiting in the queue, this system is studied as a Markov process. This type of queueing system is encountered in the operations of an automatic storage retrieval system. The steady-state probability vector is shown to be of (modified) matrix-gcometric type. Efficient algorithmic procedures for the computation of the rate matrix, steady-state probability vector, and some important system performance measures arc developed. The steady-state waiting time distribution is derived explicitly. Some numerical examples are presented.


Introduction
We consider a discrete time single-server queue in which customers arrive according to a Markovian Arrival Process (MAP).Service is provided in groups of at least one but not more than M customers at a time, where M, a fixed constant, gives the upper limit for the server capa- city.The service times of each group are assumed to follow a discrete phase type distribution whose parameters may depend on the group size.When the server becomes free, that person starts to serve the first M of the waiting customers if the number waiting is at least M; otherwise, he or she initiates a service at the next arrival for group of service size i, with probability O i, 1 <_ <_ M. With probability 1-O i, the server remains idle.Note that 0 M 1. 1This research was supported in part by Grant No. OGP0006584 from the Natural Sci- ences and Engineering Research Council of Canada to A.S. Alfa and National Science Foundation Grant No. DDM-9313283 to S. Chakravarthy.The set of probabilities 0 plays an interesting role in this type of system.The functional form of 0 is usually non-decreasing in most real life situations and the important aspect that affects the system performance is whether it is a non-decreasing convex, concave or step function.For example, a convex function reflects the behavior of a very conservative server who is not overly anxious to start service until the number of the customers waiting is very close to M. A concave function reflects the behavior of an aggressive server who is anxious to start service without keeping customers waiting for too long.In both cases, the server uses judgment to decide when to start serving groups less than M. On the other hand, a step function refers to a deterministic case where the server knows when and when not to start service.Later in the paper we investigate the impacts of these functional forms of 0i on the performance measures of the system.
This type of queueing system is encountered in an automated storage and retrieval system where the material handling system has a capacity of carrying a maximum of M items.A deci- sion has to be made whether to proceed and attend to the waiting items when there is less than M waiting.Usually the service time involves the long haul to the appropriate region in the racks plus the time per item to be placed in rack positions.Hence, the service times may be dependent on the group size.
Another application of this model is in the package delivery service by a courier company.Packages to be shipped by such companies are often time sensitive.For packages with common destination, the company would like these in batches of reasonable size.The company may thus decide, for economic reasons, to have a minimum number of packages to form a batch before they are shipped, and a maximum number at which the packages must be shipped.In between those two numbers the dispatcher has to decide when to ship the packages.In this example the vehicles are the servers and the packages are the customers.
Neuts and Nadarajan [10] considered problems of this class with many servers.They assume Poisson arrivals and exponential services with parameters independent of the group size.In their mode] a server is not allowed to serve less than a, and not more than b customers.Sim and   Templeton [11] considered the same type of model but allow a probability distribution for the batch size (between a and b).Later Chakravarthy and Alfa [1] extended that model to include a Markovian arrival process and also introduced a probability O of serving the waiting customers (a<i<b).
Earlier Medhi and Borthakur [6] considered the case involving two servers and general bulk service.Medhi [5] later extended this model to the case involving c-servers.

The Model
The arrival process is assumed to be a MAP described by two substochastic matrices D o and D of dimension n.The MAP was introduced in Neuts [9] as a generalization of the Poisson pro- cess which is well suited for matrix analytic and numerical investigations.In conjunction with the research reported by Lucantoni, Meier-Hellstern and Neuts 3], Lucantoni [2] suggested a con- venient notation which is better suited for a general discussion than the one which was originally used.
A highly accessible discussion of the MAP, with many examples, may be found in Lucantoni [2].A partly expository paper discussing how the MAP can be used qualitatively to model point processes with certain "bursty" features is given by Neuts [7].The following is a brief informal description of the discrete MAP, (DMAP), which should be adequate for the purposes of this paper.Let D, D o / D 1) be an irreducible stochastic matrix of order m and let D O and D 1 be two substochastic matrices whose sum is D and such that the matrix I-D o is nonsingular.
At any transition from state to state j of the Markov chain with transition probability matrix D, that transition is an arrival epoch with conditional probability (D1)ij[Dij] -1 and is not an arrival epoch with the complementary conditional probability (Do)ij[Dij] -1.Equivalently, we may consider the discrete-time Markov renewal process embedded at arrival epochs and with transition probabilities defined by the sequence of matrices Q(k) [D0] k-1D1, for k _> 1.
The DMAP is a discrete-time point process generated by the transition epochs of the above Markov renewal process.The simplest DMAP is the Bernoulli process for which the matrices D o and D 1 are scalars, respectively given by, 1-p and p.The point process is then the sequence of epochs of successes in independent Bernoulli trials.
The stationary vector r of the Markov chain described by D O + D 1 satisfies the equation: r(D 0 + O1) r, re 1 (1) where e is a column vector given e [1,1,...,1]', where denotes the transpose operation and rDe is the probability that, in the stationary version of the arrival process, there is an arrival at an arbitrary time point.Correspondingly, * rDie is the expected number of arrivals per unit of time and is also referred to as the fundamental rate of the MAP.
Service is provided in groups of sizes i, 1 _< _< M, with probability 0 i.The service for a of size is of phase type with representation (fli, Si) of order n i.For later use, let S be group such that Sic + S e.
Note that (#)-1 i(i_ oi)-le, where is the service rate of the system when a group of customers is served.
Let us define some notations that will be used in the sequel.The symbol (R) stands for the Kronecker product.That is, X (R) Y stands for the matrix made up of blocks XijY.For more details about Kronecker products see Marcus and Minc [4].
Suppose U and V are square matrices of order m and n respectively.
g and V, denoted by U 4 Y W, is an (n + m) (n + m) matrix with Wij Uij(i j 1,..., m) (2) (4) In a partitioned form we write Hence The variable v is defined as v n I -t-n 2 + n 3 +... + n M and ej, is a column vector of appro- priate dimension with 1 in the jth position and 0 elsewhere.

The Markov Chain
The system described in Section 2 can be studied as a Markov chain.The state space for this system can be described by A-{(j,k),0_<j_<M-l,l_<k_<n} {(i, jl,j,k, j3), i_>0, The states and their description are given in Table 1 below.

Description
The server is idle with j customers waiting and the arrival process is in phase k.
The server is busy serving a group of J2 customers with their service in phase J3, (i M) + Jl, customers are waiting and the arrival process is in phase k.
Note that the number of customers waiting is organized in groups of equal sizes M and the remaining Jl less than M.This is because the server will serve in groups of size M whenever there is at least that number waiting.
The matrix P is of the Quasi-Birth and Death type.The steady state vector x partitioned as x Ix*, x(0), x(1), x(2),...] and associated with this matrix is obtained by solving xP x and xe-1. Theorem: Proofi The matrix A A o + A + A 2 is block circulant and its invariant probability vector is simply given as --[b,b,...,] where b is the invariant probability vector of G- It is know from Neuts [8] that A2e > A0e if the matrix P is positive recurrent.Applying this condition, and after routine algebraic manipulations, the stated result follows.
Noting that the matrix P has a structure that satisfies both the GI/M/1 and the M/G/1 paradigms as described by Neuts [8], we choose the GI/M/1 paradigm to take advantage of the special structure of the matrix A 0. The vector x(i + 1) can be determined recursively as x(i / 1) x(i)R, for _> 0, where R is the minimal non-negative solution to the matrix quadratic equation: R A o + RA 1 + R2A2.
(25) (26) (27) 4. The Matrix R By partitioning R into blocks of nu nu matrices and noting that the first (M-1) blocks of A o are zeros, R can be rewritten as: Hence, R can be computed using matrices of smaller dimensions as follows Lemma 1: We have Rj-R, 1 <_ j <_ M, where R is the solution to R 1 G 0 + R1G1 + RIMG2 + R1M + 1G3" (30) Proof: Equation (28) can be rewritten as: R (G O -1-RMG2 -t-RMR1G3)(I-G1)-1 (31 where I-G 1 is nonsingular.Replacing Rj with R in (29), it is straightforward to show that Rj as given by Rj -[(G O + R1MG2 + R1M + 1G3)(I-G1)-1]j (32) satisfies Equation (28) and Equation (29).Using the fact that the irreducibility of A guarantees a unique solution to Equation (28) and Equation (29) (see Neuts [8]), the stated result follows.Equation ( 28) can be further simplified by exploiting the structure of Go, G1, G 2 and G 3. Both G O and G 1 are block diagonal matrices and G 2 and G 3 have the first (u-riM) blocks of columns of zeros.Thus R 1 can be written as 0 where Rjj is of dimension (nnj) x(nnj) and RiM is of dimension (nnj) x(nnM).The computation of R 1 is carried out as follows.First the diagonal (block) matrices Rjj, 1 _ j _ M-1 are computed.Secondly, RMM is evaluated and then RiM 1 _ j _ M-1 is evaluated.
The following lemma gives relevant expressions.
4. The probability that the number of customers waiting is less than kM-1, Yk, is given by Yk E ki= lx(i-1)e.

Numerical Examples
6.1 PH/Geo/1 System In this example we consider a PH/Geo/1 system, and our interest is to determine the effects that various theta values would have on the system performance measures.
In all cases we take O.
where > 1 is a lower threshold before service is even considered.
Example 1" In this example we take M 8, D O and D are as follows, 0.2 0.6 0.06 0.14 C 0.3 0.4 0.09 0.21 The vector of the service rate is [0.66, 0.4356, 0.2875, 0.1897, 0.1252, 0.0827, 0.0546, 0.036].The arrival rate for this is 0.2655.Figures 1 and 2 display the two performance measures, mean queue length and server idle probability as a function of the group size.An examination of these figures reveal the following observations.For high traffic where PM-0.036, the shortest queue length is for 0i(1).Also Y0 is smallest when group size is M.
Example 2: In this example we consider a system with low traffic by taking #M-0.101.
Arrivals have the same MAP representation as in Example 1 and the vector of the service rate is [0.75, 0.5625, 0.4219, 0.3164, 0.2373, 0.1778, 0.1335, 0.101].Figures 3 and 4 display the two per- formance measures, mean queue length and server idle probability, as a function of group size.An examination of these figures reveal the following observations.For high traffic where #M-- 0.101 the shortest queue length is for 0i(1).Also Y0 is smallest when group size is M.
Figures 1 through 4 show that the best performance was found when theta values are 0/1.
The results indicate that a shorter queue length results when there are specific guidelines that tell the server when to serve the customers.

MAP/PH/1 System
Since we determined from the PH/Geo/1 examples that theta values of 0/1 type result in the shortest queue values, we will only consider 0/1 type values for the thetas in the remaining exam- ples (see Appendix B).
We take0j-0, j<iand0j For example, if j-3, the server cannot start serving unless there are at least 3 customers waiting for service.Negative Binomial distribution which is a special case of a discrete phase distribution and will be indicated in the examples.The Negative Binomial distribution is the discrete analog of the Erlang distribution.As group size increases so does the average service time.Our selection of the Negative Binomial distribution for service is due to its more common occurrence in real life.This is not a limitation of our model and our model allows a general phase type distribution.

Arrival Rates
Example 3: In this example we take phase service to be of the same dimension n 2 phases, (i 1,...,10).M 10, fl= fli [1,0] (i 1,...,10), and theta values are of the form 0/1 type.Our interest in this example is to study the system for the case where phase service has the same dimension for the different group sizes even though the mean service times are increasing with group size.
Table 2 shows the service times for each group size i; p-.1 is the mean service time and (r. is the variance of the service time. The results are shown in Figure 5.The mean queue length increased as the group size increas- ed.According to this example, the smallest mean number in the queue resulted when the group size was 1.The smallest Y0 value (probability the server is idle) also resulted when the group size was 1.As the traffic intensity increased (from 0.2 to 0.9) so did the mean number in the queue and Y0 (see Figure 6).Y0 increases initially as increases and then it decreases as increases.In this example, both the mean and variance of the service time increased monotonically with i.  Example 4: In this example, we allow the dimension of the phase service to vary with the group size and take n i--l, i=1,2,3; n i=2, i--4,5,6,7; n i=3, i=8,9,10.M--:10, fli-1, i= 1,2,3; fli [0,1], i= 4,5,6,7; fli [1,0,0], i= 8,9,10, and theta values chosen as 0/1 type.
In this example, the mean of the service time for each group size is the same as that for the corresponding group size in Example 3.However, the variance of the service time is non-mono- tonic.
Table 3 shows the service times for each group size and the corresponding mean and var- iance.The results are shown in Figure 7.According to this example, the minimum mean number in the queue resulted when the group size was 1 for low to medium-high traffic.However, for high (0.9) traffic intensity, the minimum mean number in the queue resulted when the group size was 4.
Note that the service time of group size 4 has the smallest variance.This seems to have a major impact on the performance measures at a high traffic level.Increasing the mean service time did not necessarily result in an increase in mean queue length, the mean queue length either decreases or did not increase at an expected rate when the variance of the service time is decreased.This is to be expected.As shown in Figure 8, in all cases server idleness probability was smallest with a group size of 1, with the lowest value obtained with a traffic intensity of 0.9.

The Stationary Waiting Time Distribution
In this section we present the equations for obtaining the waiting time distribution.computational aspect is intensive and is a subject of future work.

The
Before obtaining the waiting time distribution we need to set up some notations.Let yj denote the probability vector that an arriving customer finds the server idle with j waiting in the queue.Similarly y. .,(i),i>0, 0<j<M-1, I<j'<M is defined as the probability vector ,3 that an arriving customer finds the server busy with j' customers and there are iM + j customers waiting.Then we have y*--,x(O)(I(R)D1) i>_O,  y(i) .x(i)(I (R) D1) >_ O, (51) and * is the arrival rate.
Further, let P(r,k) denote the matrix of dimension n n whose (i, j) th entry is given by Pij(r,k) Pr{N(k) r, Jk j N(O) o,g o-i} where N(k) number of arrivals during (O,k] and J phase of the arrival process at time .It is easy to verify that P(r,k)-P(r,k-1)D0+P(r-l,k-1)D1, r>k+l where I r-O P(0, r) O, r:/:O where ($(i),T(i)) of order nMi is the representation of the convolution of discrete PH-distribu- lions, all with the same representation ([JM, SM), and $(i)TU2-Ul l(i)T(i) E P(s,k 1 u2)Dle" r + j +s + 2, k >_ 3. s--O Proof." Follows immediately on noting the following.
1.For w(0), i.e., the customer gets into service immediately upon arrival: In this case the server is idle when the customer arrives and there are j, 0 < j < M-1 customers waiting in the queue.The server therefore starts to serve j + 1 customers with probability O j + 1" 2. For Wl(k), k >_ 1" the customer arrives to find the server idle with j, 0 < j <_ M-2 cus- tomers waiting.The customer waits, and enters service after k units such that during the first k-1 units of time r, 0 _< r < M-2-j, customers arrive, and at time k an arrival triggers the service.
3. For both w21(k), k >_ 1 and w22(k), k _> 2: the customer arrives to find the server busy with Jl customers and there are j customers waiting in the queue.
(a) For w21(k), k >_ 1: there are two mutually exclusive events to consider.

i)
The current service lasts exactly k units, and 0_< j <_ M-2.The customer enters service k units after arrival with another customer who arrives at time k.
During the first k-1 units of wait r, 0 < r < M-2 j, customers arrive.ii) Here j, 0 _< j _< M-1, customers are waiting in the queue.The current service lasts k units.During the first k units at least Mj customers arrive.
(b) For w22(k), k>_2 the server ends service at time u, l<u_<k-1.The server is idle between time u and k and initiates a service at time k with an arrival.During the interval (0, u] there are r, 0 < r < M-2-j, arrivals and during the interval (u,k-1] there are s, 0 < s _< M-2-r-j, arrivals. 4. For both w31(k), k _> 2 and w32(k), k > 3: the customer arrives to find the server busy with Jl customers and there are iM + j, >_ 1, customers waiting in the queue.
(a) For w31(k), k >_ 2" the current service lasts u units and the service time of the groups of size M lasts k-u units.The customer enters service k units after arrival at which time another customer arrives.We have to consider two mutually exclusive events.

i)
During the first k-1 units r, 0_<r<M-2-j, customers arrive, 0_<j_< M-2, ii) During the k units wait at least M-1-j customers arrive, 0 _< j _< M-1.(b) For w32(k), k >_ 3" the server ends the current service at time ul and serves the groups of M customers in u 2-ul, l_<u l_<u2-i, i+l<u2_<k-1 units of time.The customer enters service at time k after arrival with another customer who arrives at time k.During the first u 2 units r, 0 _< r _< M-2-j, customers arrive and during the next k-1u 2 units s, 0 _< s _< M-2-r-j, customers arrive.

Table 1
States and their description State

Table 2
Service times for group size i

Table 3
Service times for group size i