ON THE CALCULATION OF STEADY-STATE LOSS PROBABILITIES IN THE GI / G / 2 / 0 QUEUE

ABSTILCT This paper considers methods for calculating the steady-state loss probability in the GI/G/2/0 queue. A previous study analyzed this queue in discrete time and this led to an efficient, numerical approximation scheme for continuous-time systems. The primary aim of the present work is to provide an alternative approach by analyzing the GI/ME/2/0 queue; i.e., assuming that the service time can be represented by a matrix-exponential distribution. An efficient computational scheme based on this method is developed and some numerical examples are studied. Some comparisons are made with the discrete-time approach, and the two methods are seen to be complementary.


Introduction
In this paper we consider methods for calculating the steady-state loss probability in the GI/G/2/0 queue. In a previous study of this system [3], numerically exact probabilities were obgained for discrete-time systems in which both the interarrival time and the service time take finite integer values, and this led to an effective numerical approximation scheme for continuoustime systems.
Our primary aim in the present work is to provide an alternative method for solving continuous-time problems. We do this by using a matrix-exponential representation of the service-time distribution; i.e., by studying the GI/ME/2/0 queue. A matrix exponential density function f(x) has the form f(x) a exp(T_ x)s__, where _a is a row vector, _s is a column vector and _T is a matrix, and complex entries are allowed in a, s and T. An equivalent characterization is 1This author's research for this paper was supported by the British Council that f(x) has a rational Laplace transform. The matrix exponential distribution is identical in form to the widely-used phase-type distribution (see Nests [5]), except that in the latter case the entries in a,s and T are subject to additional restrictions. Both matrix-exponential and phasetype distributions are powerful models for approximating general probability density functions in queueing and related processes. In the next section, we give a detailed steady-state analysis of the GI/ME/2/0 queue, leading to two computational schemes for the calculation of loss probabilities. This is followed, for comparative purposes, by a brief description of the discrete-time model referred to above. Some preliminary, numerical studies are then presented to illustrate both the discrete-time and matrix-exponential approaches and to make some comparisons between them.
2. The Queueing System GI]ME[2]O We begin with some notation" A(x) distribution function of interarrival time, b(x) probability density function of service time, t n arrival time of customer n, Ptoss steady-state loss probability; the desired value, u n number of busy channels at time t n / 0, n service time of customer n if u n 1, (n1), (n 2) residual service times in channels 1 and 2 respectively, if u n n ((n 1), (n2)) is a vector.
We make the following assumptions: --2; in such cases Both the mean interarrival time a and the mean service time s are finite; moreover a > 0. b(x) takes the following matrix-exponential form: where it suffices to assume that Re(,i) > 0 and the right-hand side of (1) should be nonnegative for all x > 0.
Obviously the sequence (tn,n) is a homogeneous Markov chain. One could easily derive integral equations for its steady-state distribution, but we shall use Cox's principle of complex probabilities [4], following which we can write down steady-state equations as if i were positive and bik > 0 and summing up to 1. Thus each customer arrival originates a trial with outcome probabilities bik so that in the (i,k) case the service time is the sum of k independent exponential phases having parameter i" Therefore, instead of the variables n, one may consider the embedded Markov chain ('n, in, in) in which n has the same meaning as above. If u n 1 then, for the customer in service at time t n +0, the parameter of its exponential phase is i (i.e., having index in) while Jn is the number of residual phases including the current one. If un-2, (2) index the parameters of the exponential (1), i)), j (j(l) j2)) where i and z then n (z n phases in channels 1 and 2 respectively, while j(n 1) and j(n 2) are the corresponding numbers of residual phases in channels 1 and 2.
i(1)i,i(2)i,,Fi(1),i(2),j(1)_ j,, j(2) j,, Thus, for example, fij is the probability of completing j exponential phases within an interarrival interval, the parameter of the exponential law being h i. The above transition diagram implies the following steady-state equations: i1,'2,31, J2 m How many variables are involved? The system (4) contains (' ri) 2 variables z(2, il,i2, J1,J2) but the symmetry condition 1 z(2, i2, il, J2, Jl) z(2, il, i2, Jl, J2) For many widely used distributions, both the series in (15) converge sufficiently rapidly. For example, if the service-time distribution has an increasing failure rate, then X n and Yn tend to zero in an exponential manner. In general, one cannot assure good estimates for X n and Yn, but both certainly vanish as n goes to infinity. Thus, in either case, one needs a finite number N of iterations for a given admissible error in the value of Ploss" How many elementary operations are needed? It can be seen that, at each iteration, an (i1, i2) variable at step n is connected to (i',i2) and (i1,i') variables at step n-1. It therefore suffices to take the number L of elementary operations each consisting of a multiplication and an addition, required for each iteration, to be given by t--1 In practical computations, very large values of r can be expected to occur infrequently, otherwise the model of b(x) would have a low chance of being assured statistically, so that an estimate of 2NL would seem to be quite satisfactory. Of course, it remains to be seen how large N, the required number of iterations, will be in practical applications.
Mixture of Exponentials: In this case, b(x) takes the form m i=1 We can simplify the notation as follows: x>0.
Each iteration requires L 2m 3 elementary operations, not counting the symmetry reduction.
"Ityper-Erlangian" Model: In this case b(x) takes the following form" Equations (13) Each iteration requires L-2r 3 elementary operations, not counting the symmetry reduction.
It is shown in [3] that the Ck (k >_ 1) satisfy the following system of equations: The quality

Numerical Examples
In this section we describe some preliminary numerical investigations which allow us to illustrate, and make some comparisons between, the matrix-exponential approach and the discrete-time approach to the GI/G/2/0 queue. The examples involve the following family of distributions B n (n 1,2,...) having the density: bn(x on(1-cos(27rx))exp(-x), 0 x _ n (29) where c,, (1 + 4)(1 exp( n)) 1.
For convenience, as n goes to infinity we replace Bn, b n and c n by B, b and c respectively.
The densities corresponding to finite n are truncated versions of b(x). The latter, which is shown in Figure 1, is given by Asmussen and Bladt [1] as an example of a density function that does not belong to the class of phase-type distributions but does belong to the larger class of matrixexponential distributions. The matrix-exponential representation is: b(x) a exp(T_ x)s_ where 2ri-1 0 0 -c/2 _a - (1,1,1), _T-0 -2ri-1 0 _s -c/2 0 0 -1 c and, in the form of (21), b(x) is a mixture of exponentials in which m-3, 1-1-2ri, 2 1 + 2ri, 3 1 and b 1 -c/(2)tl) b 2 -c/(22) b 3 c/. 3. Bn/B/2/0 and Bn/Bn/2/0 have been studied. The first of these queues has Ploss values given by the well-known Erlang Loss formula; this queue is of interest here partly to validate the computational method and also to study the rate at which the algorithm converges. Obviously, one expects the rate of convergence to depend upon the traffic intensity s/a, and this can be seen in Figure 2 where the estimate ptoNss! of the steady-state loss probability is plotted against the iteration number N; i.e., modifying (15), Next we consider the queue Bn/B/2/O for a range of n values, in each case calculating loss probabilities by both the matrix-exponential model and the discrete-time model. In the latter case, the distribution B14 (i.e. corresponding to L-14) was taken to be a sufficiently accurate representation of B, the truncated tail probability having a value c that is less than 10 -6. The number of time points rn used in the approximation took values up to a maximum of 7,000, and the algorithm was terminated when the value of Ploss calculated for the approximating problem had been determined to an accuracy of 0.0001. The results for the two different approaches, using programs running on a Vax 11/780 mainframe computer, are compared in Table 1. For the Bn/Bn/2/O queue, in which both the interarrival time and the service-time distributions have the form B n with n finite, the matrix-exponential method is no longer applicable as an exact method. In this case, only the discrete-time approach has been used and the results are given in Table 2.
The results in Table 1 show good agreement between the two methods. The Ploss values for the matrix-exponential method are, of course, exact to the prescribed accuracy, while those for the discrete-time method are approximate but with errors that are quite low. For the examples using 7,000 time points, the errors are about 0.3%, which is comparable with results for other systems described in [3]. In general, the discrete-time method was heavy in its use of computing time (see, for example, the c.p.u, times listed in Table 2) while, for the matrix-exponential method, in each case the c.p.u, time was negligible. It is interesting to note that in the latter case the number N* of iterations could be predicted quite closely using the graph in Figure 3 for the M/B/2/0 queue. Of course, the distribution B is defined precisely and is of low dimension (r 3). In a more general case, it might be necessary to fit an appropriate ME density to represent an actual service-time distribution, for example by using the approach described in [2]. In this case, the speed of the matrix-exponential method would clearly depend upon the dimension of the ME representation used, and the accuracy would depend upon how well the actual service-time distribution was represented by the fitted ME density. In the absence of an appropriate ME representation, the discrete-time method is sufficiently flexible to give acceptable approximations, as shown in Table 2 for the Bn/Bn/2/0 queue.
Finally, it should be noted that the form of the discrete-time distribution fitted to the B n densities in this paper, is not necessarily the most appropriate one for all applications. In the trade-off between computing time and accuracy, simpler discrete-time approximations could be devised by combining fewer time points, and hence shorter computing times, with a less precise fit to the distribution 's shape, for example by matching moments. The implications of this trade-off will be the subject of further research.

Concluding Remarks
In this paper we have presented a computationally efficient method for calculating steadystate loss probabilities in the queue GI/G/2/0 when the service-time distribution can be represented by a matrix-exponential distribution. We have also considered an alternative discrete-time approximation method. Preliminary numerical studies suggest that the two approaches are complementary in that the preferred choice in any application will depend largely upon the nature of the service-time distribution and how accurately it can be modelled by each of the two approaches.