MULTI-THRESHOLD CONTROL OF THE BMAP / SM / 1 / K QUEUE WITH GROUP SERVICES

We consider a finite capacity queue in which arrivals occur according to a batch Markovian arrival process (BMAP ). The customers are served in groups of varying sizes. The services are governed by a controlled semi-Markovian process according to a multithreshold strategy. We perform the steady-state analysis of this model by computing (a) the queue length distributions at departure and arbitrary epochs, (b) the LaplaceStieltjes transform of the sojourn time distribution of an admitted customer, and (c) some selected system performance measures. An optimization problem of interest is presented and some numerical examples are illustrated.


Introduction
One of the most popular and effective ways to model the flows of messages in modern communication networks or jobs in a production and manufacturing process is to use the batch Markovian arrival process (BM AP ).The BM AP is a rich class of point processes that includes many well-known processes such as Poisson, Switched Poisson, Markov-modulated Poisson, and P H-renewal processes.In recent years, there has been a constant and growing interest in the investigation of queues with BM AP input [1,6,17,41].
The behavior of such queues is described usually in terms of continuous or discretetime Markov chains.These fall under the so-called M/G/1 type Markov chains [44] or the quasi-Toeplitz Markov chains [23].These chains have two important properties.The first one is the Toeplitz-like structure for the (block) transition probability matrix governing the queueing system under study.That is, the (i, l) th block matrix P i,l of the one-step transition probability matrix of the denumerable Markov chain corresponding to the transition from state i to state l depends only on l−i, but not on i and l separately for i > i 0 where i 0 ≥ 0 is some integer.The second property is the so-called non-skip free to the left property.That is, P i,l = 0, for 0 ≤ l < i − 1.
The presence of these two properties enables one to study the stationary distribution of the chain using matrix analytic methods.We refer the reader to [1,6,17,23,28,32,34,36,37,44] for latest developments in this area.
The "non-skip free to the left" property is an important one in the development of the theory of Markov chains of the M/G/1 paradigm.If we do not assume this property (e.g., queues where negative arrivals take place or where services are in groups) the investigation of the Markov chain governing the corresponding queueing system becomes far more difficult (see, for example, [2,3,4,35]).However, there are at least two interesting cases where the analytical investigation of the chain is still possible when the skip free to the left property is not satisfied.
The first case arises in the investigation of queues with disasters.Disaster is a special case of negative arrivals of customers.A disaster causes the removal of all customers (including those in service) from the system instantaneously.For more information on such queueing systems we refer the reader to the survey paper [2].A Markov chain that describes the queue with disasters has non-zero blocks P i,0 for any i ≥ 0. Hence, such chains do not belong to the class of quasi-Toeplitz type Markov chains.However, the stationary distribution of the queue length can be obtained in a rather nice analytical form (see, for example, [27,29,33]).
The second case is where the customers are offered services in groups of varying sizes.In the context of a finite capacity queueing system with finite buffer of size K, a different type of service scheme in which services are offered to groups of varying size, ranging from a predetermined threshold L to the maximum buffer size K, was introduced in [7].The pre-assigned number L ≥ 1, called the threshold, operates as follows.An idle server finding fewer than L customers in the queue remains idle.However, when i, L ≤ i ≤ K, customers are present, the idle server initiates a service for the entire group.If a service has to be initiated through an arrival of a batch (for example, i jobs, i < L, are in the queue with the server being idle, and the arriving batch has at least L − i jobs) only a group of size no larger than K enter into service, and any remaining customers in the batch are considered lost.Service schemes of this type in the context of finite capacity GI/P H/1 and M AP/G/1 with single arrivals, and BM AP/M/c models were investigated in the papers [7,8,9,10,11,12,13,14,15,16].Some potential applications of this type of service mechanism in manufacturing processes were outlined in those papers.
Optimization problems in queueing theory play an important role in practice.They can be broadly divided into two groups.One group deals with static optimization, where some system parameters must be tuned to provide the desired quality of service to customers.For example, A.K. Erlang investigated the M/M/N/0 model and when the arrival and service intensities are fixed, he calculated the minimal number of channels needed to guarantee that the loss probability does not exceed a prespecified value.
The other group of optimization problems that falls under dynamic optimization involves the control of some parameters of the queueing system such as the arrival rate, the service rate, and the number of channels so as to optimize a given objective function.This class of problems is more complicated and challenging as the objective function is a (nonlinear) function of the steady-state probabilities and for most interesting practical queueing models these probabilities will not be known explicitly.The traditional way of solving such optimal control problems in queues is to use Markov decision processes.However, this approach is not very powerful even in the case of classical queueing systems.Furthermore, in the case of Markov chains describing queues with BM AP input and SM -services there is the curse of dimensionality problem.Thus, it appears that when parametric strategies of control such as threshold, hysteretic, and randomization are used in queueing models, the problem of finding an optimal control is better solved by a direct approach.This approach involves an iterative process consisting of three stages: (1) Computation of the stationary queue length distribution for a given set of values for controlled parameters, (2) Evaluation of the objective function, and (3) A search for optimal values for the controlled parameters using an efficient heuristic approach.
Such an approach was used in analyzing the BM AP/G/1 system with a multi-threshold service rate control with N service modes in [26], for the BM AP/P H/1/K queue with hysteretic control in [5], for the BM AP/G/1 queue with hysteretic control in [30], and for the BM AP/SM/1/K queue with hysteretic control in [31].
In this paper, we extend this approach to the BM AP/SM/1/K system with group services.This generalizes the results of the paper [5] in two aspects.The first one is the assumption of availability of J, J ≥ 2, service modes (in [5] J = 2).Secondly, we assume SM -type service time distributions while in [5] P H-type services were used.Furthermore, here we make a more realistic assumption about the service switching mechanism in that we allow the switching only at the beginning of a service epoch.In [5] the service rates can be changed at any epoch using a hysteretic type control.
The paper is organized as follows.In Section 2 the mathematical model and the service control mechanism are described.The stationary queue length distribution embedded at service completion epochs is investigated in Section 3 and the stationary queue length distribution at an arbitrary time point is obtained in Section 4. In Section 5 some key system performance measures describing the queueing model are presented along with their formulas.The stationary sojourn time distribution in the system of an admitted customer is derived in Section 6. Section 7 contains a brief description of an optimization problem and some illustrative numerical examples are presented in Section 8.

The Mathematical Model
We consider a single server queue with a limited buffer of size K. Arrivals of customers occur according to a batch Markovian arrival process (BM AP ), a tractable class of Markov renewal process.Customers in a batch are admitted to the extent of buffer availability and all others are considered lost.That is, we allow "partial admission" of a batch.We assume that the behavior of the BM AP is governed by the directing process {ν t , t ≥ 0}.This process is an irreducible continuous-time Markov chain with the finite state space {1, . . ., W }. Suppose that the matrix D 0 governs the transitions corresponding to no arrivals, and the matrices D k , k ≥ 1 govern the transitions corresponding to arrivals of batch size k, k ≥ 1.By assuming D 0 to be a nonsingular matrix, the interarrival times will be finite with probability one and the arrival process does not terminate.Thus, we see that D 0 is a stable matrix.Let N t denote the number of arrivals in (0, t] of the {ν t } process.In the sequel, we need the matrices P (n, t), n ≥ 0, t ≥ 0 where the (j, k) th entry of the matrix P (n, t) is defined as where For use in the sequel, let 0 and 1 denote, respectively, the row and column vectors of 0's and 1's with appropriate size.The matrix D(1) is the infinitesimal generator of the chain {ν t , t ≥ 0}.Let u denote the stationary vector of this generator.That is, uD(1) = 0, u1 = 1.The (group) arrival rate, λ g , and the fundamental rate (or the average intensity), λ, of the BM AP are defined as: and λ = uD (1)1.
For full details on BM AP we refer the reader to [32,42,44] and for a review and recent work on BM AP we refer the reader to [17].
We assume that the service is offered to groups of varying size i, L ≤ i ≤ K.Here L is some pre-assigned number, L ≥ 1.If at the completion of a service fewer than L customers are present, the server waits until the queue length reaches at least L and then initiates a service to all those customers present in the system.
In order to attract and serve customers efficiently it may be necessary to control the service rate (or the service time distribution) according to the number of customers in the queue.We assume that the (group) service times of the customers follow an irreducible finite state semi-Markov process {m t , t ≥ 0} with state space {1, 2, . . ., M}.The sojourn times of the services are given by the entries of the kernel B(t) = ||B m,m (t)|| of dimension M .In this paper we categorize the services into J, J ≥ 2, modes of operation.During the j th mode, the service times are given by the entries of the kernel m,m (t)|| of dimension M .In other words, while the services of all groups are directed by the common process {m t , t ≥ 0}, the transitions of this process are governed by different kernels depending on the type of strategy control used.We further assume that the semi-Markov services times have finite second moments.Note that independent and identically distributed service time is a special case of the SM -service.In this case, the service time of the group is characterized by the distribution function We consider a multi-threshold strategy in this paper.This is defined by the integervalued thresholds I 1 , . . ., If the number i of customers in the system at a service completion satisfies the inequality I j−1 < i ≤ I j , then the entire group of i customers is served in the j th mode, 1 ≤ j ≤ J.If the number of customers i is less than L, no service is offered until the queue builds up to L or more.In the latter case, we consider two possible variants for selecting the service mode for the next group.In Variant 1 the service for the next group will always be in mode 1 and in Variant 2 the service mode for the next group is determined at the beginning of the service according to the multi-threshold strategy.Note that in the case of MAP (which corresponds to single arrivals) these two variants coincide.Other types of variants for choosing the service mode can easily be incorporated without any additional complexity and the details are omitted.In the case when the thresholds are equal, the number of service modes is reduced appropriately.Thus, if r thresholds are equal then the number of service modes will be reduced by r − 1 to J − r + 1 and the service modes will be relabelled as 1, 2,..., J − r + 1.
In situations where there are various costs associated with the waiting of customers as well as for providing faster modes of services, one of the most popular classes of the parametric strategies is the class of the multi-threshold strategies.The optimality of such multi-threshold strategies in the class of all homogeneous Markovian strategies is proven only in some particular cases (see, e.g., [18,46]).However, such strategies are practical as well as reasonable for numerical implementations.Various queueing systems, such as M/M/N [19], GI|M |1 [22,38], M |G|1 [20,21,39], BM AP |G|1 [26], BM AP |SM |1 (with retrials) [25], and BM AP |SM |1|N [24], with multi-threshold strategies of control have been analyzed in the literature.For the type of group services introduced in [7], to our knowledge this is the first time such a multi-threshold policy is considered.

Embedded Queue Length Distribution
Let t n denote the epoch of the n th service completion; i n the number of customers in the queue at t n +0, i n ≥ 0; ν n the state of the BM AP process, {ν t }, at t n , 1 ≤ ν n ≤ W ; and m n be the state of the semi-Markovian process, {m t }, that governs the service process at t n + 0, for 1 ≤ m n ≤ M, n ≥ 1.When the server becomes idle immediately after a service completion due to not having enough customers in the queue, we assume that the phase at t n + 0 will be the phase of the next service to be initiated.That is, the phase of the service process will be frozen at the instant when the server becomes idle.We can modify this scheme to allow the phase to be chosen according to some initial probability vector, but this will not be addressed in this paper.
Let the thresholds I 1 , . . ., I J−1 be fixed and predetermined.It is easy to verify that the process {(i n , ν n , m n ), n ≥ 1} is a three-dimensional Markov chain.The assumptions of our model imply that this Markov chain has a stationary state distribution.

Define
π(i, ν, m) = lim Enumerating the states of the Markov chain (i n , ν n , m n ) in lexicographic order, we form the vectors π(i), of dimension W M, for 0 ≤ i ≤ K, of stationary probabilities.
We now define a number of auxiliary quantities that we need in the sequel.Let Ω (j) i denote the matrix of probabilities that during the service of a group of customers in the j th service mode exactly i customers arrive.That is, and Let X (i) denote the f irst passage probabilities of going from level 0 to level i or higher.That is, the components of X (i) give the probabilities that the process {N t , ν t , m t } reaches level i or higher for the first time from level 0. It is easy to verify that X (i) is given by to be the matrices of probabilities that the process {N t , ν t , m t } reaches level l starting from level m.These matrices are given by and Lemma 1: The vectors π(i), i ≥ 0 satisfy the following equations: for Variant 1 (where service is always in mode 1 once the group size builds to L or more); and for Variant 2 (where the mode of service is determined at the beginning of a service after the group size builds to L or more) the equations are ) The quantity χ(l) is equal to Proof: The proof follows directly from the law of total probability.Theorem 1: The stationary probability vectors π(i) are calculated as follows: where the vector v = (v 1 , . . ., v J ) is the unique solution to the system v = vA, (3.12)

.13)
The matrix A is defined as follows.In the case of Variant 1, the entries A rj are given by: with For Variant 2, the entries A rj are of the form:

16)
Proof: We will prove for Variant 2 as the proof is similar for the other case.Defining and substituting (3.17) into (3.8) and (3.9) we get expressions (3.10) and (3.11).Thus, to complete the proof we have to show that the vectors v j , 1 ≤ j ≤ J satisfy system of equations given in (3.12) and (3.13) for Variant 2. To this end, we first multiply equation (3.10) by and add over 0 When this is added to the sum of equations ( 11) over I j−1 ≤ i ≤ I j , we get where the matrices A rj are defined by formula (3.16).The stated result follows immediately. Remarks: (1) Note that the vectors π(i), 0 ≤ i ≤ K can be calculated directly by solving the system of linear equations in (3.6) and (3.7) (or (3.8) and (3.9)).However, when J < K, use of Theorem 1 will reduce the computational efforts.
(2) The entries of the matrices A rj of dimension W M have a very nice probabilistic interpretation.Partitioning the matrix A rj into W blocks of M by M matrices, the (l, l ) th entry of the (k, k ) th block matrix gives the transition probability that the next service will start in mode j with phase l and at that instant the arrival process will be in phase k given that the current service is in mode r with phase l and the arrival process is in phase k.

Stationary Queue Length Distribution at Arbitrary Time Points
Let z(0, i), 0 ≤ i ≤ L−1, be the row vector of size W M whose entries z(0, i, j, k) give the steady-state probabilities that at an arbitrary time there are i customers in the queue with the arrival process in phase j and the server became idle from service phase k.
be the vector whose entries z(n, i, j, k) define the steady-state probability that at an arbitrary time there are i customers in the queue with the arrival process in phase j and the current service for a batch of n customers is in phase k.To derive an expression relating the steady state probabilities at an arbitrary time and at embedded epochs, first we need the following lemma dealing with the mean interdeparture time.Lemma 2: Suppose that τ denotes the mean interdeparture time.It is calculated as where Proof: The proof follows immediately from the law of total probability on noting that the entries of the matrix Ĝ(L−i) give the (conditional) mean time until the beginning of the next service for a group of L or more customers given that that the previous service left behind i, 0 ≤ i ≤ L − 1, customers in the queue.
Theorem 2: The stationary-state probability vectors z(n, i) in Variant 1 are obtained as follows: where and B (t) is the diagonal matrix with the diagonal entries defined by the vector B (j) (t)1, 1 ≤ j ≤ J, and the indicator function χ(n) is as defined in Lemma 1.
In Variant 2, the formulas are corrected by means of replacing ∇ (1)

5).
Proof: Proof follows from a classical argument based on the key renewal theorem (see e.g.[40]).
Remark: Note that the arbitrary time stationary probabilities depend implicitly on the type of variant considered for initiating a service while waiting for the queue size to build up to L or more, through the vectors π(l).Corollary 1: Let η i be the vector such that the (j, k) th component gives the probability that at an arbitrary time there are i customers in the system and that the arrival and service processes are respectively in phases j and k.Then we have Note that z(0, i) = 0 for i ≥ L and in the above summation when the lower bound is greater than the upper bound the value is set to zero.Corollary 2: The probability vector, q(n), n ≥ L, of seeing n customers at the beginning of a service is obtained as (4.7)

Selected System Performance Measures
In this section we give a number of system performance measures and their respective formulas useful in qualitative interpretation of the model.
(1) The probability that an arriving customer will be lost is given by (2) Let ζ j , for j ≥ 0, denote the probability that exactly j customers are lost at an arrival epoch.Then we have where λ g is the (group) arrival rate of the BM AP .
(3) The mean number of lost customers at an arrival epoch, µ NL , is calculated as The throughput of the system is given by λ(1 − P reject ).
(5) Denoting by θ 0 the fraction of time the server is idle and by θ r , 1 ≤ r ≤ J, the fraction of time the server is busy serving customers in r th mode, we have ) ) in the case of Variant 1 and in the case of Variant 2.
(6) The mean number of customers in the queue, µ QL , is given by (5.8) (7) Suppose that S rj , for r = j, 1 ≤ r, j ≤ J, denotes the average number of service switches per unit of time.Then from the probabilistic interpretation of the quantities v r and A rj , we have where A rj is as defined in (3.14) for Variant 1 and (3.16) for Variant 2.

Sojourn Time Distribution
Let V(x) be the vector distribution function of the sojourn time in the system of an admitted customer at an arrival epoch.We will call this admitted customer a tagged customer.Partitioning V(x) into V(x) = (v 1,1 (x), ..., v 1,M (x), ..., v W,1 (x), ..., v W,M (x)), the entry of v j1,j2 (x) gives the distribution function of the sojourn time of the tagged customer given that the arrival process is in state j 1 and the current service is in phase j 2 .Recall that if the server is idle then the new service will start in phase j 2 .Let e −sx dV(x), Re s > 0, be the Laplace-Stieltjes transform of V(x).
Before we derive an expression for v(s) for the tagged customer, we need the following result.Suppose ξ l,r , 0 ≤ l ≤ L − 1, L ≤ r ≤ K, denotes the sojourn time from the epoch when the server is idle with l customers waiting in the system to the epoch when the service starts for a group of r customers.

Proof:
The proof follows immediately on noting that In order to derive an expression for v(s) for the tagged customer, we need to consider two cases.
Case 1: Suppose that the tagged customer is part of a batch of k customers at the epoch when the server is idle and when i, 0 ≤ i ≤ L − 1, customers are waiting in the queue.In this case we need to consider the following three scenarios depending on the size of k: In this case the tagged customer and the others enter service immediately and the service time distribution function is given by B (χ(k+i)) (t).
Case 1B: Suppose that k is such that 1 ≤ i + k ≤ L − 1.In this case the tagged customer waits for a random duration, ξ k+i,r until the number of customers in the system reaches r, L ≤ r ≤ K, before entering into service.The service time of the tagged customer will then be given by B (χ(r)) (t).
Case 1C: Suppose that k is such that i + k > K.In this case the tagged customer along with K −i−1 customers will enter into service immediately and their service time is governed by B (J) (t).
Case 2: Suppose that the tagged customer arrives as part of a batch of k customers at the epoch when the server is busy with m customers and that there are i customers waiting in the queue.In this case, we need to consider the following two scenarios.
Case 2A: Suppose that k is such that i + k ≤ K. Suppose that l customers arrive during the residual service time of the current service.Note that the current service is governed by B (χ(m)) (t).If i + k + l ≥ L, then the service time of the tagged customer is given by B (χ(i+k+l)) (t).If i + k + l < L, then the tagged customer waits for a duration ξ i+k+l,r , L ≤ r ≤ K before entering service which is governed by B (χ(r)) (t).
Case 2B: Suppose that k is such that i + k > K.In this case the tagged customer along with K −i−1 customers will enter into service immediately and their service time is governed by B (J) (t).
Combining all of the above scenarios, the following relation holds good for x > 0: Defining the following theorem is easily verified.

Theorem 3:
The vector Laplace-Stieltjes transform v(s) of the sojourn time of an admitted customer in the case of Variant 2 is given by where Note: The expression in (6.4) can be used to calculate the mean sojourn time of an admitted customer in the case of Variant 2.
Corollary 3: The (marginal) Laplace-Stieltjes transform v(s) of the tagged customer is given by v(s) = v(s)1.(6.5) Remark: Theorem 3 can easily be modified for Variant 1 and the details are omitted.

An Optimization Problem
In this section we consider an optimization problem.Let c j , 0 ≤ j ≤ J, denote the cost per unit of time of service in mode j.Note that when j=0, the server is considered to be idle.Let d 1 denote the holding cost per customer per unit of time of waiting in the queue and d 2 the cost per lost customer per unit of time.Then the optimization problem of interest is given by min I1,...,IJ−1 where P reject , θ j , and µ QL are as given in (5.1),(5.4)-(5.8).
Finding an optimal solution in the set of all multi-threshold policies is very complex.Furthermore, the solution of this problem is complicated due to the fact that the objective function is known only implicitly in terms of the steady state measures.Hence, developing a very clever and efficient heuristic algorithm and the numerical implementation of this algorithm require a very careful and detailed analysis.In the next section we provide some interesting numerical examples.

Numerical Examples
In this section we will discuss some illustrating numerical examples.We consider three arrival processes with the following BM AP representations {D k }.For all these BM AP s it is easy to verify that λ g = 2.8261 and λ = 4.2391.While the first two arrival processes correspond to renewal processes, the third arrival process corresponds to a correlated process with the correlation between the successive interarrival epochs given by 0.12795.The standard deviation of these three arrival processes are, respectively, 0.25021, 0.53731, and 0.43652.

BMAP with Positive Correlation:
We take J = 3.That is, the system has at most three operation modes.The service times are of semi-Markov type.For all three operation modes, the transition matrix P of the embedded Markov chain for the directing process {m t , t ≥ 0} and the semi-Markovian kernels B (j) (t), 1 ≤ j ≤ 3, are taken to be where the distribution functions B In the following we consider Variant 1 for the service mechanism.That is, the service initiated through an arrival of a customer will always be in mode 1.
In Tables 1 through 3 that follow we list the values of θ j , 0 ≤ j ≤ 3, P = P reject , and µ QL for various combinations of the thresholds I 1 and I 2 , for the three arrival processes.Note that I 3 = K = 10.Also, when I 1 = I 2 we have only two service modes, namely mode 1 and mode 3, and when I 1 = I 2 = I 3 , there is only one service mode which is mode 1.From these tables we notice the following observations.
• For all I 2 and for the three arrival processes considered, the mean queue length (µ QL ), the loss probability (P reject ), and the fraction of time the server is busy in service mode 1 (θ 1 ) appear to increase as I 1 increases.The rate of increase decreases as I 1 increases.It is interesting to note that as I 1 increases, the mean queue length for Erlang arrivals appears to dominate the other two arrival processes for all values of I 2 .
• For all I 2 and for the three arrival processes considered, the fraction of time the server is idle (θ 0 ), the fraction of time the server is busy serving in mode 2 (θ 2 ), and the fraction of time the server is busy in service mode 3 (θ 3 ) appear to decrease as I 1 increases.The rate of decrease appears to decrease as I 1 increases.
• With respect to the measure, θ 2 , there appears to be a cut-off point, say, I * 2 (I 1 ), such that for I 1 ≤ I 2 < I * 2 , Erlang has the largest value for this measure and hyperexponential has the least value.For I 2 ≥ I * 2 , hyperexponential has the largest value and Erlang has the least.For example, when I 1 = 2, we have I * 2 = 2 and when I 1 = 4, I * 2 = 6.The same type of phenomenon appears to hold true for the measure θ 3 .
• While for all values of I 1 and I 2 , the measure θ 2 appears to decrease with increasing variance of the arrival times, the measures θ 0 and P reject appear to increase with increasing variance of the arrival times.However, more experimentation is needed to see how correlation plays a role.
• We see that the optimal value 9.085900 for the cost criterion in (41) occurs at I 1 = 10 and I 2 = 10 for Erlang arrivals; the optimal value of 9.621708 occurs at I 1 = 4 and I 2 = 9 for the hyperexponential case; and for BM AP with positive correlation the optimal value of 9.36168 occurs at I 1 = 6 and I 2 = 7.It is interesting to note that for this choice of parameters, there is no control mechanism required for the Erlang arrivals; however, for the other two arrival processes, control strategy yields a better solution.This indicates that whenever the arrival processes have a larger variation or are correlated the situation will be quite different and needs a careful analysis.

Table 1 :
Cost function and its components for various values of I1 and I2 for Erlang arrivals

Table 2 :
Cost function and its components for various values of I1 and I2 for hyperexponential arrivals

Table 3 :
Cost function and its components for various values of I1 and I2 for positively correlated arrivals