APPROXIMATIONS TO PERFORMANCE MEASURES IN QUEUING SYSTEMS

Approximations to various performance measures in queuing systems have received considerable attention because these measures have wide applicability. In this paper we propose two methods to approximate the queuing characteristics of a GI/M/1 system. The first method is non-parametric in nature, using only the first three moments of the arrival distribution. The second method treads the known path of approximating the arrival distribution – by a mixture of two exponential distributions – by matching the first three moments. Numerical examples and optimal analysis of performance measures of GI/M/1 queues are provided to illustrate the efficacy of the methods, and are compared with benchmark approximations.


OPSOMMING
Benaderings tot verskeie prestasiemaatstawwe in toustaanstelsels ontvang beduidende aandag weens die wye toepasbaarheid daarvan.In hierdie artikel word twee metodes voorgestel om die toustaaneienskappe van die GI/M/1-stelsel te benader.Die eerste metode is nie-parametries van aard en gebruik slegs die eerste drie momente van die aankomsverdeling.Die tweede metode volg die bekende roete om die aankomsverdeling te benader deur die eerste drie momente te pas, deur middel van 'n kombinasie van twee eksponensiёle verdelings.Die doelteffendheid van die metodes word aan die hand van numeriese voorbeelde en optimale ontleding van die prestasiemaatstawwe van GI/M/1stelsels bewys, en word vergelyk met benaderings.a a a a a a a a a a a a a a a a a a a

INTRODUCTION
The growth of queuing theoretic applications has been phenomenal, ranging from communication and multimedia systems to inventory and reliability theory.This has led to a sustained interest in methods to evaluate performance measures in queuing theory.In the case of non-Markovian queues, the computations of these measures involve the arrival and/or service distributions explicitly.However, in practical applications like management, optical, and communication networks, the specific forms of these distributions might not be known.At best, one might only be in possession of the moments of the underlying distribution.There are a number of cases where the moments of a distribution are easily obtained, but theoretical distributions are not available in closed forms [9].Alternatively, from the sample data observed, efficient estimators for the various moments of the underlying distributions can be calculated.Thus, computation of performance measures based on the first two or three moments of the arrival and/or service distributions is very useful.Whitt [16], in a classic exposition, discussed approximations using extremal distributions giving the upper and lower bounds for the performance measures in a GI/M/1 system.Smith [13] proposed a two-moment approximation for the probability distribution of M/G/1/K systems, and extended it to the analysis of M/G/1/K queuing networks.Sohn & Lee [14] conducted a Monte Carlo simulation in order to study the relation between various performance measures in a G/G/1 queue.Recent work on such systems with working vacations for the server have immense applications in ATM machines and internet systems, such as optical nets, electric nets, and communication nets [8,4,2].In these applications, the arrival epochs could be observed or, at worst, simulated.Our motivation in this paper has thus been to obtain approximations to the performance measures of a GI/M/1 system using only the first three moments of the arrival distribution, without explicit recourse to the arrival distribution We will discuss our problem with specific reference to a GI/M/1 queuing system, even though our methods also work in a similar way for other non-Markovian queues.Consider a GI/M/1 queue whose traffic intensity is ρ= E (service time)/E (arrival time), L is the expected equilibrium queue length, and σ is the steady state probability that a customer will have to wait to begin his service.It is well known that where σ is the unique root in the open interval (0, 1) of the equation with µ=1/E(service time) and Ф(s), the Laplace-Stieltjes Transform of the inter-arrival distribution function, say F, given by: We note that the evaluation of the performance measures σ and L require prior knowledge of the inter-arrival distribution function F, and not just the moments of F. As mentioned in the beginning of this section, many queuing applications are likely to produce only the moments of F and not the distribution itself.Thus the problem is to find σ and L on the basis of the first few moments of F only.Whitt [16] showed that there is a considerable reduction in the range of possible values of σ and L when the third moment is also used, compared with using just two moments of F. We propose simple and accurate methods to evaluate σ and L based on the first three moments of the distribution function F in the absence of any knowledge of the form of F. In section 2, we propose a non-parametric method based only on the first three moments of F -without recourse to approximate F -by another distribution function.Numerical illustrations are provided to compare the values of σ and L, using the present method, with their exact values.The method provides exact results for certain important arrival distributions like Erlang of order 2, Coxian (K 2 ), a mixture of two exponentials, and exponential distribution.We also provide two optimisation illustrations to obtain economic performance measures in the application of GI/M/1 queuing systems.
Approximations of probability distributions by phase type distributions (by matching moments up to a certain order) have attracted the attention of researchers because of necessity, and because of their wide applicability.Pioneering work on phase type distributions and their various applications was done by Neuts [11].Among the various members of the family of phase type distributions that have been studied, mixtures of two exponential distributions (known as H 2 distributions) play a key role in many approximations used in queuing theory.These distributions are log convex in nature, and can approximate highly-skewed distributions quite accurately.In section 3, we suggest a simple nonlinear programming method in which the first two moments are matched exactly, while the third moment is matched as closely as possible.This method works in the entire region of possible values (of m 1 , m 2 , m 3 ), and provides exact three-moment match wherever possible.
The approximated H 2 distribution with the given three moments of the inter-arrival distribution is then used to calculate σ and L. Numerical illustrations are provided to validate the approximation and computation of the performance measures.

A NON-PARAMETRIC METHOD
We observe from equation ( 2) that the computation of the performance measures σ and L in the GI/M/1 system requires the use of Ф(s), the Laplace Transform of the density function f corresponding to the distribution function F. However, without prior knowledge of F, and armed only with the first three moments of F, an approximation to Ф(s) is obtained from the proposition that follows.

Proposition
Suppose that the first three raw moments m 1 (≠0), m 2 , and m 3 (3m 2 2 ≠2m 1 m 3 ) of the distribution function F exist and are known.Then the following approximation to the Laplace Transform of the distribution function F holds.

Proof
In the classical renewal theory, the renewal density m(t) of a renewal process with interval density f(t) satisfies the integral equation: Applying Laplace Transform to both sides of (6), we obtain the Laplace Transform of m(t) as: where Ф (s) is the Laplace Transform of the density function f.Now Ф (0) =1 and   �1 − ()�| =0 =  ́(0) =  1 is non zero.This means that the denominator 1-Ф (s) of ( 7) has a simple zero at s=0.Thus we may approximate m * (s) by a rational function of the form where A, B, and s 0 are constants determined as follows.Assuming the existence of moments of the density function f(t), we can express Ф (s) as where m 0 =1 and m n is the n th order moment about the origin of f.Using ( 9) in ( 7) and ( 8), and comparing the coefficient of s 0 , s 1 and s 2 on both sides, we obtain (using some algebra) the values of A, B, and s 0 as given in ( 5).This completes the proof.
Note1: For the approximation to hold, it is necessary that s 0 ≤0.Simple calculations show that this condition implies that Ф 2 >2 and (see Figure 1), where C 2 is the squared coefficient of variation of inter-arrival time, Ф 2 = C 2 +1, and

Figure1: Feasible regions for the approximation
Note2: It is worth mentioning that the error in the approximation ( 4) is of o(1) as s →0.
Note3: The condition that s 0 is non-positive, which is necessary for the approximation to hold, is satisfied by many standard arrival distributions: uniform, gamma, mixed exponential, lognormal, Coxian (K 2 ), mixture of Erlangian, and Weibull.Also s 0 is nonpositive for truncated normal and inverse Gaussian probability density functions under certain conditions.Using (4) in equations ( 2) and ( 1) with the values of m 1 , m 2 , and m 3 , the performance measures σ and L were obtained immediately.To illustrate the efficiency of the proposed method, we present (in Table 1) the values of σ and L computed using (4) for certain choices of the set {m 1 , m 2 , m 3 }.In order to compare the approximations with exact values, we have considered the values of the moments of commonly used arrival distributionsnamely, gamma (Table 1) and PH4 distributions (Table 2).The values of σ and L are calculated for various values of traffic intensity ρ.Eckberg [5] specified upper and lower bound distributions that yield the maximum and minimum mean queue length in a steady state among all inter-arrival time distributions with two and three moments specified.Using these distributions, Whitt [16] calculated the maximum relative error for σ and L using the formula where L l and L u are the minimum and maximum values of L using the lower and upper bound distributions.Table 1 also presents the upper and lower bounds for σ and L given two and three moments, and the corresponding maximum relative errors specified in (10).It can be seen that our method captures the values of σ and L with low relative errors.

INTER-ARRIVAL DISTRIBUTION APPROXIMATION BY MATCHING MOMENTS
Approximating general distributions by phase type distributions is important in queuing theory because their structure leads to Markovian state description and consequently analytical tractability.Although several phase type distributions have been used in the literature, two distributions that are prime candidates for such approximations (because of their simplicity and suitability) are mixtures of two exponentials (H 2 ) and Coxian (K 2 ) distributions.On this point, we use the former for analysis, as these distributions provide a fairly accurate match when the C 2 is large -which is true for arrival distribution in a queuing system.Furthermore, in some of the examples discussed by Whitt [17], the maximum relative error (MRE) when two moments are fitted was found to be 200 percent, while working with mixtures of exponentials reduced the MRE to 50 percent.Specifying the third moment reduced the MRE to 5 percent.Thus, a three-moment match using H 2 distributions for inter-arrival distributions seems to provide useful results.
Using an empirical study, Bere [3] showed that when the service time distribution is approximated using the first two of its moments, the third moment has a considerable effect on the probability distribution of the number of customers in an M/G/1 queue if . He also showed that the probability distribution of the number of customers and the average number of customers in λ(n)/G/1/N system are highly sensitive to the third moment of the service time distribution if C 2 >1.Altiok [1], in justifying the inclusion of the third moment in matching, refers to Bere's empirical work.Whitt [15] empirically showed that the effect of the third moment on the average number in the system in a GI/G/1 queue becomes considerable as C 2 increases.If the service time distribution has C 2 <1, the impact of the third moment is not significant [3,15].Since the present work deals with matching an H 2 distribution with three moments, we confine our attention to the range C 2 >1 only.
It is well known [15] that three numbers .Thus, in the region where exact three-moment matching is not possible, researchers have used adhoc methods to find the approximate H 2 distribution.These regions are clearly shown in Figure 2. The parameters p, μ 1 , and μ 2 of H 2 (see ( 11) below) when (m 1 , m 2 , m 3 ) falls in Region-I, and where the three moments can be matched exactly, are well documented.However, when these moments fall in Region-II, such that the three moments cannot be matched exactly, methods suggested in the literature are recipes in nature.Lopez-Herrero [10], in the absence of information on service distribution, used the maximum entropy principle approach to estimate the true distribution of the number of customers served during the busy period in an M/G/1 retrial system.Whitt [15] suggests that "if m 3 turns out to be too small when attempting an H 2 -fit, one procedure is to replace m 3 by something slightly larger than (3/2) m 2 2 /m 1 ".Altiok [1] suggests the use of m 3 = 3m 2 2 /2m 1 + εm 1 3 but does not indicate how to calculate the perturbation parameter ε.In the following algorithm, we propose a goal programming procedure of matching the first two moments exactly, and matching the third moment as closely as possible in Region-II.This procedure subsumes Region-I as a particular case.We note that Region-III is an infeasible region for the existence of (m 1 , m 2 , m 3 ).

Algorithm
Given the first three moments of a distribution function F (say, m 1 , m 2 , and m 3 , which are assumed to be finite), the parameters p (0<p<1), μ 1 , and μ 2 of the H 2 distribution given by (11) -whose first three moments either match m 1 , m 2, and m 3 exactly, or match the first two moments exactly, and match the third moment as closely as possible in the sense of squared differences -are given by the following algorithm: Step 1: Find optimal p using the Golden section method [12,6] to solve the following one dimension optimisation problem: Step 2: Using the value of p found in step 1, compute: The steps of the algorithm are justified using the following arguments.Consider the following probability density function of a mixture of two exponentials (H 2 distribution): In the following, we do not consider the trivial cases of p=0, 1 as they lead us to exponential density.
The first three moments of the density function above are: From (12) we have Substituting ( 15) in ( 13), and after simple algebra we have two cases.First we consider We know that Also, μ 1 >0 implies that  17) in (15) results in The third parameter p is obtained by matching the third moment as closely as possible.Thus, we minimise ( ) subject to: In the second case, we set Using algebra similar to the first case results in and the third parameter p is determined by solving It is easily seen that both cases lead to the same result, but with the roles of p and 1-p interchanged.
In Table 2 we continue with PH4 inter-arrival distribution.To illustrate this, we fit H 2 distribution for this distribution by matching the moments specified by the algorithm.The performance measures σ and L obtained by using the fitted H 2 distribution in steps 1 and 2 are given in Table 2. Attention is drawn to the MRE for these measures when three moments are matched [16] and to the actual relative error in using the approximated H 2 distributions.Also note that these two methods improve in accuracy for increasing ρ as expected.

Illustration 1
In practice, queuing managers are generally interested in optimising the model parameters under their control by minimising the operating cost or maximising the business profit.In the first illustration, we will be interested in obtaining the optimal service rate in a cost minimisation problem for a GI/M/1 queuing system.The objective cost function consists of two components: the cost due to customers waiting in line (known as the delay cost), and the service cost rate.Thus the cost function to be minimised is given by: where λ and µ are the arrival and service rates respectively, W is the expected waiting time of a customer in the system, c 1 is the expected cost per unit time of a customer's wait, and c 2 is the service cost rate.Using Little's formula, (24) reduces to The optimal µ* of the objective function above was computed using our moments matching method (introduced in section 3) by assuming the first three moments of the arrival distribution only, and the cost rates.However, in order to compare our results with the exact values, the moments were chosen to correspond to Coxian (K 2 ) and Inverse Gaussian distributions commonly used in queuing theory.The results are presented in Tables 3 and 4.
When the Coxian arrival distribution was used, our method provided the exact values of µ*.
In the case of Inverse Gaussian distribution, the relative errors were significantly small.Table 3.The optimal service rate µ* with Coxian inter-arrival distribution.(The optimal µ* computed, using our method and using the Coxian distribution, match exactly)

Illustration 2
In the second illustration, we consider an optimisation problem in a GI/M/1 queue with working vacation for the server.These problems have wide application in Internet systems such as optical, electrical, and communication nets [8].We consider a single server queuing system that has the general arrival process.The working vacation and vacation interruption are connected, and the server enters into vacation when there are no customers, and it can take service at the lower rate during the vacation period.If there are customers in the system at the instant of a service completion during the vacation period, the server will π return to the normal working level regardless of whether the vacation ends.Otherwise, it continues the vacation.The performance measure L, the mean queue length, and P(J=0) and P(J=1) -which are the state probabilities of a server in the steady state -have been derived by Li et al. [8].We refer the reader to their paper for the relevant expressions.Li et al. [8] considered the problem of optimising the service rate η during the server's vacation period for a given cost structure.Let c w represent the unit time cost of every waiting customer, and c 1 and c 2 the service costs per unit time during the normal working level and vacation period respectively.The expected net cost function to be optimised can be seen to be min: Z = c w L + c 1 µ P(J=1) + c 2 η P(J=0) (26 where µ is the service rate during the service period.The optimal service rate η* was computed using our non-parametric method of section 2 for certain values of the model parameters and cost parameters.We have used the Coxian arrival distribution and its moments (as used in Illustration 1) to obtain the optimal service rate η*.However, the Inverse Gaussian distribution could not be used, as the objective function (26) loses its convexity and becomes monotonic.We have used Erlangian of order 2 (used by Li et al. [8]) in its place.In Figures 3 and 4, we present the values of η versus the associated cost.The optimal η* and the corresponding cost obtained using our method are very close to the exact values.This study introduces two methods for evaluating performance measures in a GI/M/1 queuing system in the absence of information on the arrival distribution, and when only the first three moments are known.The first method is non-parametric as it does not use the distribution function, whereas the second method uses an H 2 distribution obtained by moment matching procedure.This procedure involves the computationally economical Golden section method.Note that a Coxian (K 2 ) distribution is also a good phase type distribution to consider.It is worth pursuing the regions (Ф2, Ф3) in which each of these approximations scores higher than the others in terms of relative errors.The usefulness of the methods in optimisation procedures has been illustrated with examples.

Figure 2 :
Figure 2: Region of three-moment matching (Region-I: Exact three-moment match possible; Region-II: Exact three-moment match not possible; Region-III: Infeasible region for three moments of a distribution function)