A queueing system with queue length dependent service times, with applications to cell discarding in ATM networks

A queueing system (M/G1,G2/1/K) is considered in which the service time of a customer entering service depends on whether the queue length, N(t), is above or below a threshold L. The arrival process is Poisson, and the general service times S1 and S2 depend on whether the queue length at the time service is initiated is <L or ≥L, respectively. Balance equations are given for the stationary probabilities of the Markov process (N(t),X(t)), where X(t) is the remaining service time of the customer currently in service. Exact solutions for the stationary probabilities are constructed for both infinite and finite capacity systems. Asymptotic approximations of the solutions are given, which yield simple formulas for performance measures such as loss rates and tail probabilities. The numerical accuracy of the asymptotic results is tested.


Introduction
1.1 Background Queueing systems arise in a wide variety of applications such as computer systems and communications networks.A queueing system is a mathematical model to characterize the system, in which the arrivals and the service of customers (users, 1This research was supported in part DOE Grant DE-FG02-96ER25168.D.I. Choi was supported by the Korea Research Foundation.Current address of D.I. Choi: Department of Mathematics, Halla Institute of Technology, San 66 HeungUp- Li HeungUp-Myon, WonJu-Shi, KangWon-Do, 220-712, South Korea.36 D.I. CHOI, C. KNESSL, and C. TIER packets or cells) occur randomly.The customers arrive at the facility and wait in the queue (or buffer) if the server is not available.If there are many customers in the queue, they may suffer long delays which causes poor system performance.Thus, the arrival rate or the service rate may need to be controlled to reduce delays.These sys- tems may be represented by queueing systems with queue-length-dependent arrival rates or service times.That is, if the queue length exceeds a threshold value, the arri- val rate may be reduced (e.g., overload control), or the service rate may be increased (e.g., the cell discarding scheme of Choi and Choi [2]).Many schemes of traffic con- trol in ATM (asynchronous transfer mode) networks have been analyzed using such threshold-based queueing systems (see Choi and Choi [2], Gong, et al. [5], Li [8], Sriram, et al. [9-11], and Takagi [12]).
In this paper, we analyze a queueing system with queue-length-dependent service times.Customers arrive at the queue by a Poisson process, and there is only one ser- ver.The service times of customers depend on the queue length.Concretely, we specify a threshold value L for the queue.If the queue length at a customer service initiation is less than the threshold L (respectively, greater than or equal to the threshold L), the service time of the customer follows a distribution with probability density function bl(.) (respectively, b2(.)).We believe that our analysis can be extended ot the case of multiple thresholds.Both infinite (M/G1,G2/1) and finite capacity (M/G1,G2/1/K) queues are considered.
The analysis of this queueing system was directly motivated by the cell-discard- ing scheme for voice packets in ATM networks (see Choi and Choi [2], and Sriram, et al. [9, 11]).In Sriran, et al. [11], a system with deterministic service times was pro- posed, in which voice packets are divided into high and low priority ATM cells.The cells arrive as a concatenated pair (i.e., two cells per arrival), and the threshold L and the capacity K are measured in terms of cell pairs.The cell discarding occurs at the output of the queue and immediately prior to transmission, based upon 'the total number of cell pairs in the queue.If this number is less then L, then the next pair of cells is transmitted (no cells are discarded).However, if the queue length is greater than or equal to L, then only the high priority cell is transmitted.Thus, the low priority cell is discarded.The analysis in Sriram, et al. [11] is based on numerically solving an embedded chain formulation of the problem.The system is studied only at service time completions.Our results, when specialized to G D (deterministic service), are directly applicable to this model.
We analyze this queue system using the supplementary variable method.We first consider the case of an infinite capacity queue, and obtain an explicit formula for the steady-state queue length distribution Pn" When the service times have different exponential distributions, the queue length distribution has a simple, closed form.
We also compute asymptotic approximations to the queue length distribution Pn for various choices of the system parameters.Next, we examine the finite capacity queue, and again obtain explicit expressions for Pn, and in particular the probability PK that the queue is full and the probability P0 that it is empty.Also, we investigate asymptotic approximations for the queue length distribution for large values of the threshold L and the queue capacity K.We show that this queueing system has very different tail behavior (and hence loss probability) than other M/G/1 type models, and that the service tails can sometime determine the tail of the queue length.
There has been some previous analytical work on queueing systems with queue- length-dependent service times (see Abolnikov, et al. [1], Dshalalow [3], Fakinos [4], Harris [6], and Ivnitskiy [7]).Harris [6] considered the MX/G/1 queue with queue- length-dependent service times.In particular, if there are customers in the queue, the service time of the customer starting service has a general distribution depending on i.By using the embedded Markov chain method, Harris [6] derived the probabili- ty generating function for the queue length at the departure epochs.However, the ob- tained probability generating functions contains infinitely many unknown constants.
A closed form was obtained only for some special service times of two types.Fakinos [4] analyzed the G/G/1 queue in which the service discipline is last-come first-served and the service time depends on the queue length at the arrival epoch of each custom- er.Abolnikov, Dshalalow and Dukhovny [1] considered queues with bulk arrivals (i.e., compound Poisson input) and state-dependent arrivals and service.Assuming that the state-dependence applies only when the queue length is below a critical level, the authors characterize the generating function (for both transient and steady-state cases) of the queue length probabilities in terms of the roots of a certain equation.Ivnitskiy [7] also considers a model with bulk, state-dependent arrivals and state-de- pendent service.Using the supplementary variable method, he obtains a recursion relation for the Laplace transforms of the transient queue length probabilities.For a very good recent survey of work on state-dependent queues, we refer the reader to Dshalalow [3].
The model here is a special case of that studied in Ivnitskiy [7].However, we are able to give more explicit analytical expressions, from which we can easily obtain asymptotic expansions for tail probabilities and loss rates.These clearly show the qualitative dependence of these performance measures on the arrival rate and the ser- vice distribution(s).In particular, we show that the tail behavior of the model with threshold L is much different than the tail behavior of the standard M/G/1 and M/G/1/K models.

Problem Statement
We let N(t) be the queue length at time t, including the customer in service, and let X(t) be the remaining service time of the customer currently in service.Customers arrive according to a Poisson process with arrival rate /, and are served on a first- come-first-serve basis.There is a single server and a queue with finite capacity K.
An arrival that would cause N(t) to exceed K is lost, without effecting future arri- vals.The service time of each customer depends on the queue length at the time that customer's service begins.If the queue length at service initiation is less than L, the service time of that customer is $1, while if the queue length is greater than or equal to L, the service time is S 2. The service times S 1 and S 2 have density functions bl(.and b2(.), with means m 1 and M1, respectively.We define Pl rnl and P2 "MI" The process (N(t),X(t))is Markov, and X(t) is referred to as a supple- mentary variable.
An important local balance result can be obtained by integrating the balance equations with respect to x from 0 to cx, which leads to: Pn -t-1 (0) " J Pn(X)dx' n >_ 1.

(1.11)
In the following sections, we construct exact solutions to the infinite capacity model and then the finite capacity model.As we will show, the solution of the infin- ite capacity model can be used to construct the solution of the finite capacity model.Then we obtain simple formulas for the performance measures by constructing asymp- totic approximations to the exact solutions.

M/M1,M:/1 Queue
To illustrate the important characteristics of the solution, we first consider the case in which the service times, S 1 and 5'2, are exponential with probability density functions: The solution of Equations (1.4)-(1.8)can be constructed in a straightforward manner.
0 Theorem 1: M/M1, M2/1 Queue.Let bi(x)-]1ie-tix and Pi-A/]1i for i- 1,2, and assume that the stability condition P2 < 1 is satisfied.The marginal probabi- lities are given by" Pn POP'S, n O, 1,...,L-1, (2.11 Pn popL1-1 A + 1 2 p where -p l-pj" An interesting aspect of the result is the tail behavior as n-.The tail pro- bability has a different form depending on the parameters: We now consider the system in which the service times S 1 and S 2 have general distri- butions.For n < L, we solve Equations (1.4)-(1.6)with L-cx by introducing the generating function where Gx(z,x OG(z,x)/Ox..17) x The unknown function G(z, 0) is found by setting x 0 in Equation (2.17) to obtain: G(z, O) "PZ(Z -^1 )bl(,-,z), (2.18) Z-bl(-Az where bl(S is the Laplace transform of bl(t)..19)We invert G(z,x) to find an integral representation of the stationary probabilities in the form "Po / Z --1 x e'\(z -l)(t-X)b 1 (t)dt dz. (2.20)

P() -( z[z 1( z)]
The contour C is a loop in the complex z-plane which encircles the origin but excludes any other poles of the integrand.By setting n-1 and x-0 in Equation (2.20), we find that Equation (1.4) is satisfied.Based on the above calculation, we see that Equation (2.20)is in fact a solution of Equations (1.5)-(1.6)for L < Thus, we assume that for n < L, the solution is of the form" for some constant k, which will be determined later.

C
The quantities of interest are the marginal probabilities Pn, which are obtained by removing the dependence on x.We first compute H(z)-/ H(z,x)dx-E Pn zn-L" 0 n_L The following identity is useful in calculating the marginal probabilities: e A(x U)bj(u)du dx A 0 x By integratingH(z,x) and simplifying, we find that" (2.31) In addition, we assume that L _> 3, which is needed for the above simplification and to avoid a degenerate problem.We invert Equation (2.31) to obtain the marginal probabilities for n-L, L + 1,.. k [[ bl("-"z)-z w-1 dwdz Pn (27rit2 z] ] z b"2(,, Az) w L-lzn-L + l [w_b(A_Aw)](z_w) The contours for the double complex integral are such that [w > [z on the w-con- tour, and both are small loops about the origin.
The marginals for n < L are obtained by integrating Equation (2.21) with respect to x, which leads to" k f-l+bl(A-Az) dz, (2.33) C and, using Equation (1.4), we find that" k /(z-1)bl(,-z) l_sdz k.To complete the solution, we find the constant k using the normalization condi- tion from Equation (1.10).First, we compute" where on the contour C + we choose [w > 1.If 1 < 1, the only poles inside C + are atw-0andw-1.Since -bl(0)-m land -b(0)-Ml, We find that: where on the contour C+ we have wl >1 and wl <1+6.In 1< w] <1+6, the function w-bl(A-Aw) is assumed non-zero.
For n < L, we must compute Po + nL-11pn" We simplify Equation (2.33) as follows: Pn which leads to k / z-l _dz Combining the two sums, we find that: We simplify Equation (2.37) by noting that if Pl < 1, then: c+ By shifting the integration contour, we can remove the restriction that Pl < 1, and obtain the alternate form for P0 in Equation (2.43).
Also, the

Asymptotic Approximations
Simple formulas can be obtained by asymptotically expanding the exact integral re- presentations for the stationary probabilities given in Theorem 2. Important limits are L >> 1 and n-L >> 1. ^The asymptoti[ expansion in this limit^d epends on the location of the zeros of w-bl(-w), w-b2(-w), and poles of bl(W).Specifical- ly, we need to locate the singularities of the integrand which are closest to the origin.
Let A be the non-zero solution of: + A / eAZbl(z)dz.
We assume that bi(O are analytic for some 0 with (0)< 0, which insures the existence of a unique solution to Equation (2.50).Clearly, A > -A and satisfies > < A0ifpll.We let B be the solution of: + B / Z%(z)dz, 1) 0 which satisfies B > 0 when P2 < 1 and B0 as p2T1.The Laplace transform b1(0) may have singularities in the half-plane %(0)< 0. We assume the singularity with the largest real part is at 0--C, and that: D 0 -C (2.52) bl(0) (0 -f-C) for some constants u, D > 0 For example, if b I is exponential, then b"l(O Pl + 0' so that D 1, C #1, and u 1 in Equation (2.52).
To derive the asymptotic formula for n > L, we again view the double contour integral as an iterated integral and approximate the w integration first.The pole of w-1 that is closest to the origin is at w-1 + A/A, where A is the root of W bl (, ,w Equation (2.50).Thus, if L >> 1, We define w-1 1

54)
If n-L ((1), there is no further simplification.However, if n-L >> 1, we can re- place the integral in Equation (2.55) by an asymptotic approximation.The inte- grand is analytic at z 1 + A/A, but has a pole at z 1 + B/A > 1 and (possibly) an algebraic singularity at z 1 + C/A > 1 (see Figure 1).The asymptotic approxima- tion depends on the relative sizes of C and B.

z-plane I+B/) I+C/)
Figure 1: A sketch of the z-plane indicating the location of the poles of the integrand.
where F' is shown in Figure 3.
to obtain" -n+L-1 dw z-plane Using the substitution w (C + ,) in the integral, we find that it can be expressed in terms of the Gamma function, so that for C < B, We note that if b(O) has a simple pole at 0--C, then u-1 and the algebraic factor (n L + 1)u-1 disappears.
Theorem 3: Asymptotic approximations for M/G1,G2/1 queue.For L >> 1, and A, B, and C defined by Equations (2.50)-(2.52), the stationary probabilities have the following asymptotic approximations: For n >> l and n < L, and For L >> 1 and n-L fixed (just above the threshold), the asymptotic result is given by Equation (2.55) with k-P0 given above.We specialize the above formulas to the exponential service case with bi(t #it -tit for 1,2.For this case, we find that rn 1 1/#1, m 2 2/#12, M 1 1/#2, and M 2 2/#22.We explicitly solve Equations (2.50) and (2.51) to obtain: In addition, from Equation (2.52), we find that" C-#I D-#I v-l.
Finally, we examine the case when C , B, more precisely B-C O(L-1).For L >> 1, n >> 1 and n < L, we have the approximations from Equation (2.65).We let C B + r/L, and note that the integral in Equation (2.55) now has singularities (at 1 + B/A and 1 + C/A) that are close to each other.We consider the integral: We make the change of variables z-(1 +)(1 +), and use the approximations Here the contour P' is shown in Figure 3, and it encircles both of the singularities of the integrand.The final approximation to Equation (2.55) is, for n-L---.cx, The main quantities of interest are again the marginal probabilities.The marginals for Pn, 1 <_ n < L, and L _< n < g are given by Equations (2.40) and (2.41)   in Theorem 2, respectively.Again the constant k-P0 must be recalculated.The marginal probability PK must be computed using" PK-/ PK(X) dx" 0 We use the identities and to obtain" This result can be further simplified by noting that: The final result is: l--P1 / 1 The normalization constant P0 is determined by Equation (1.10), which we write as: L-1/pj) K-1/pj\ PK (3.6) 1-1po-+ j= I -P-"b j= L -O) -PO Substituting Equations (2.40), (2.41) (with k replaced by P0), and (3.5)into Equa- tion (3.6), we find after simplifying that" The results are summarized below.

Asymptotic Approximations
We compute asymptotic approximations for the finite capacity system as K, L-oc, and for various values of Pl and P2" The asymptotic expansions of Pn, except for P0 and PK, are the same as in the infinite-capacity systems and are given in Theorem 3 <C.
for the two cases B > Next, we expand P0 given by Equation (3.12).The asymptotic approximation for the integrals depends on the location of the singularities of the integrand that are closest to the origin, and follows closely the results for the infinite capacity system in the previous section (see Equation (2.59)).For the single integral, these poles are at w=l and w=I+A/A, where A>0 if Pl <1" The dominant poles of the double >0if integral are located at z 1, z 1 + B/A, and z 1 + C/A, where C > 0 and B < P2 X 1.Thus, for K, L >> 1, K-L > 1, and B C, we find that" uniformly in Pl and P2" We can simplify Equation (3.13) for different values of the parameters.For example, since C > 0, the last term in Equation (3.13) is always negligible and may be dropped.The simplifications of Equation (3.13) are summarized below" 1.
The asymptotic expansion of PK as K,L--c is computed by using Equation (3.5) for PK/PO and the expansions of the integrals which were derived for the infin- ite capacity case in Section 2. The results are summarized below: Theorem 5: Asymptotic expansions for the M/G1,G2/1/K queue.For K,L >> 1 and A,B and C defined by Equations (2.50)-(2.52), the stationary probabilities Pn have the asymptotic expansion from Equation (2.65) for n < i and Equations (2.63)- (2.64) for n-L >> 1.The asymptotic expansion for PK is: PK ( A )L I A(p2 --1) where the asymptotic expansions of Po are given by items 1-9 above for various values of pl and P2" 4. Discussion and Numerical Results We now demonstrate the usefulness of our results for a finite capacity system, namely the M/D1,D2/1/K queue. are: For this model, the density functions of the service times bl(X ((xml) b2(x ((x-M1).The exact solution is given in Theorem 4 with Equation (4.1) in place of the general density functions.We illustrate how to numerically compute the exact solu- tion.This calculation is only feasible for moderate values of L and K.For large values of L and K, we constructed the asymptotic approximations in Theorem 5. We demonstrate the accuracy of our asymptotic results by comparison with the exact solution.As we will see below, our asymptotic results are quite useful for moderate values of n, L and K, and are extremely accurate when n, L and K are large.
To evaluate the exact solution for Pn in Theorem 4, we must compute the com- plex integrals in Equations (3.9)-(3.12).The simplest approach is to evaluate the integrals by using the method of residues.The residue at z-0 in Equation (3.9) is difficult to compute if n is large.This is a key motivation for the development of asymptotic expansions.For moderate values of n, we compute the residues using the symbolic computation program Maple.
The double complex integrals in Equations (3.10)-(3.12)can be evaluated by rewriting the term 1/(zw) (which couples the two integrals) as:

1
-1 1 -1 (zJ. z-w-w 1-z/w-w 3-0 Using this result in the integrand in Equation (3.10) allows the double integral to be separated into: nL( / )(1/--1 j o 1 z j + L n 1...dz -7 wj + L" c c where can be identified from Equation (3.10).The sum truncates at j n-L, since for j > n-L, the z integral vanishes.FinaIly, we must evaluate both integrals by computing the residues at z 0 and w 0 for each value of j.Again the calculation is only feasible for L and n-L which are moderate in size.As before, we use Maple to perform the calculation.
The calculation of the asymptotic approximation requires computing the con- stants A and B by solving Equations (2.50) and (2.51), respectively.For determinis- tic service times, there is no singularity C; thus, C oc > B. Given the constants A and B, we evaluate the Equations (2.65) and (2.63) for Pn if n < K, and Equation (3.22) for PK" The constant P0 is given by Equations (3.14)-(3.21),depending on the values of Pl and P2" In Table 1, we consider a queue with the threshold L 5 and capacity K 10.
The arrival rate is A-1 and the service times are rn I 1/2 and M 1 1/4, so that Pl-1/2 and P2-1/4.Thus, when the queue length exceeds the threshold, the service time of the jobs entering service is half of the original service time.For these values of Pi, the asymptotic values of P0 are given by Equation (3.14).As we see from Table 1, the asymptotic expansion is quite accurate for n-0, n-2 through 4, and n>8.It is not accurate for n-land n-5 through 8.This is consistent with our results in Theorem 5, since the asymptotic result in Equation (2.65) is valid for n 1 and Equation (2.63) is valid for n-L 1.It is remarkable that our results are this accurate since n and n-L are only moderate in size.We cannot expect the asymptotic solution to be accurate for n-L, since it assumes that n-L--oc.We have chosen values for L and K that are quite small.In reality, we would expect them to be of the order 102 (see Sriram,et al. [11]).For such large values of L and K, calculation of the exact solution would prove difficult while the evaluation of the asymptotic solution is straightforward.
In Table 2, we consider the same queue as in Table 1, but now K-15.The values for A and B are the same as for Table 1 since these constants are independent of L and K.We see that the exact and the asymptotic results are numerically close to those in the previous example for n < 10, as expected.For n >_ 10, the relative error starts at about 4% and rapidly decreases to under 1%.
The example in Table 3 is a queue with threshold L-8 and capacity K-15.
The arrival rate is -1, and the service times are m I -2 and M 1 -1/2, so that Pl 2 and P2 1/2.In this example, Pl > 1, so that in the absence of the threshold, the queue length distribution would be peaked near the capacity K. Now A -0.7968 and B 2.512.The asymptotic results are quite accurate when n < L. However, when n _> L the results are not accurate.
In Table 4, we retain L-8 and increase K to 25.Now the asymptotic results are accurate to within 5% for 16 _< n_< 25.In Tables 3 and 4, the results are not accurate for L 8 <_ n <_ 14.Apparently, n-L 6 is too small a value for Equation (2.63) to be useful for these particular parameter values.In Table 5, we consider the same case as in Table 3, except that now M 1 -2/3 and hence P2 2/3.Now we have A ,, -0.7968 and B 1.144.The results are to within 7% if n-L>_4 (i.e., 12_<n_<15).In Table 6, we increase K to 25.The error is at most 7% for all n-L _> 4 (i.e., 12 _< n _< 25).Our results should be very accurate for K > 25, but it then becomes difficult to evaluate the exact solution.These numerical comparisons show that the asymptotic approximations are quite robust, and are accurate even for relatively small values of L and K. Tables 4 and 6 show that when n-25-K, the two results agree to five decimal places.Loss rates can thus be calculated to a very high precision using our asymptotic formulas.
the contours are small loops inside the unit circle.

Figure 2 :
Figure 2: The contour F.This leads to the following approximation for Equation (2.55) as n-