New Approach for Finding Basic Performance Measures of Single Server Queue

. Consider the single server queue in which the system capacity is infinite and the customers are served on a first come, first served basis. Suppose the probability density function 𝑓(𝑡) and the cumulative distribution function 𝐹(𝑡) of the interarrival time are such that the rate 𝑓(𝑡)/[1 − 𝐹(𝑡)] tends to a constant as 𝑡 → ∞ , and the rate computed from the distribution of the service time tends to another constant. When the queue is in a stationary state, we derive a set of equations for the probabilities of the queue length and the states of the arrival and service processes. Solving the equations, we obtain approximate results for the stationary probabilities which can be used to obtain the stationary queue length distribution and waiting time distribution of a customer who arrives when the queue is in the stationary state.


Introduction
The GI/G/1 queue has been given a lot of attention due to its applications to a variety of manufacturing systems, communication systems, computer systems, and other related areas.The steady-state waiting time distribution and queue length distribution are important measures of performance in this model.
Steady-state waiting time distribution in the continuoustime GI/G/1 queues can be found by numerical approximations based on the theory of Fredholm integral equations [10].
In Kim and Chaudhry [11], a discrete-time version of the distributional Little's law is established.Based on this law, the queue length distribution for the discrete-time GI/G/1 queue may be obtained from its waiting-time distribution.
Neuts [12] approximated the general distributions of service time and arrival time with phase-type distributions and solved the steady-state distribution using the matrixgeometric approach.
In this paper, we consider a single server queue in which the service time and interarrival time are assumed to have constant asymptotic rates.Abbreviating "constant asymptotic rate" as "CAR, " we may denote the resulting queue as a CAR/ CAR/1 queue and any distribution with a constant asymptotic rate as a CAR distribution.In practice many distributions such as exponential, Erlang, hyperexponential, and gamma satisfy this requirement.However the heavy-tailed arrival time distributions which are often used in queuing analysis are not CAR.In order to be still able to apply the results in this paper, we may attempt to approximate the non-CAR distributions by CAR distributions.For example, in the case when there is a value  * such that the probability of the given random variable  larger than  * is negligible, we may 2 Journal of Probability and Statistics approximate the distribution of  by a CAR distribution by changing the rate of the distribution of  at time  for all  >  * to be all equal to the rate at time  * .
To analyze the CAR/CAR/1 queue, we may first discretize it by segmenting the time axis into a sequence of equal intervals of length Δ.We next derive a set of equations for the stationary probabilities of the queue length and the states of the arrival and service processes.By solving the equations, we obtain approximately the stationary probabilities.Each of the probabilities may be found more accurately by an extrapolation of the probability on the values of Δ.
From the stationary queue length distribution, we have also obtained the waiting time distribution and sojourn time distribution of a customer who arrives when the queue is in a stationary state.
The approach which is closely related to the methods in this paper is the polynomial factorization used in [7].Instead of taking note of the three variables given by queue length and the states of the arrival and service processes at the end of the time interval of length Δ, Haßlinger inspected the queue length  and the time  until the next arrival or next departure, only at time slots when an arrival occurs or a service is completed.Thus Haßlinger only needed to use a 2-dimensional state representation instead of the 3-dimensional state representation used in this paper.An important assumption used by Haßlinger is that the stationary probability of (, ) is given by     where  and  are parameters, and   is another parameter which depends on .
A major contribution of this paper is the introduction of an easily implementable algorithm for solving the equations for the stationary probabilities of the queue length and the states of the arrival and service processes.This algorithm may be used to investigate other queuing systems which cannot be represented as continuous-time Markov chains.For example, Yang et al. [13] studied a model given by an unreliable M/M/1 queue with a multistate server whose service rate deteriorates due to the shocks which occur randomly with random magnitudes.This algorithm has been successfully used to investigate the model given in Yang et al. [13] with the distribution of the interarrival time or service time changed to a fairly general distribution (see [14]).

Derivation of Equations for the Stationary Probabilities
In this section, a set of equations for the stationary probabilities of the queue length and the states of the arrival and service processes in the discretized CAR/CAR/1 queue is derived.First let () be the probability density function (pdf) of the service time and   the interval (( − 1)Δ, Δ] for  = 0, 1, 2, . ... Furthermore let where  is large enough such that Suppose a service starts at time  = 0.Then, for sufficiently small Δ, the probability that the service will be completed in the interval  1 is approximately  1 Δ, and, given that the service is not completed in  1 ,  2 , . . .,  −1 , the probability that the service will be completed in   is approximately   Δ,  = 2, 3, 4, . . .where   =   for  ≥ .
For the arrival process, let () be the pdf of the arrival time.Furthermore let where  is large enough such that Suppose a customer has arrived at time  = 0. Then the next customer will arrive in the interval  1 with an approximate probability  1 Δ, and, given that the next customer does not arrive in the interval  1 ,  2 , . . .,  −1 , the probability that the next customer will arrive in   will be approximately   Δ for  = 2, 3, 4, . . .where   =   for  ≥ .
Given that a service starts at a time in  0 , we may define the state number   of the service process at the end of   as = 0; or (ii) the service ends in   , for  ≥ 1; or (iii) the server is idle in   , min (, ) , if the service does not end in   ,  ≥ 1.
Next, given that a customer arrives at a time in  0 , we may define the state number   of the arrival process at the end of   as Let   be the queue length at the end of   and h  = (  ,   ,   ).We may refer to h  as the vector of characteristics of the queue at the end of   .
Let  ()   be the probability that, at the end of   , the number of customers in the system is  (including the customer that is being served), the service process is in state , and the arrival process is in state , where  ≥ 0,  ∈ {0, 1, 2, . . ., } and  ∈ {0, 1, 2, . . ., }.Assume that (7) exists.To find the   , we first make the following observations.
Suppose at the end of  −1 the system is not empty, and the service and arrival processes are in state  − 1 and  − 1 at the end of  −1 , respectively.Then only one of the following events can occur in the next time interval   : (a) a customer enters the system with the arrival rate   * , and, at the end of   , the vector of characteristics becomes h  = ( + 1,  * , 0); (b) a customer leaves the system with the departure rate   * , and h  = ( − 1, 0,  * ); (c) no customers enter or leave the system, and h  = (,  * ,  * ).
Here,  * = min(, ) and  * = min(, ).However if the system is empty at the end of  −1 , and the arrival process is in state −1, then either one of the following events may occur in   ; (d) a customer enters the system with arrival rate   * , and h  = (1, 0, 0); (e) no customers enter the system, and h  = (0, 0,  * ).
By setting  −1 = 1,  − 1 =  − 1 = 0 and letting event (b) occurs in   , we get When  → ∞, (8) yields, In general, for a given value of h  , we can likewise find the combinations of h −1 and the event in   which lead to h  , and obtain an equation similar to (9).Thus we can obtain the following equations: When the queue length is  = 1, When the queue length is  = 2, the expressions for  0 , 1 ≤  ≤ , are the same as those given by ( 13), (14), and (15).Other   when the queue length is  = 2 can be computed from the equations below: When the queue length is  ≥ 3, the values of all the   (except  0 ) can be computed using ( 21) to ( 24) and ( 13) to (15), whereas those of  0 can be computed using the following equations:

Stationary Queue Length Distribution
Before solving (9) to (27) in Section 2, to obtain the stationary queue length distribution, we may first let   ,   ,   ,   and   be constants and introduce the following notations: With the above notations, ( 12) to (18) in the case when  = 1 can be represented as and ( 19) to (24) together with ( 13) to (15) in the case when  = 2 may be represented as Furthermore ( 25) to (27) together with (21) to ( 24) and ( 13) to (15) in the case when  ≥ 3 may be represented as It can be shown that, from the set of equations given by (30), we can get By substituting the expression of  1 given by ( 33) into (31), and solving for  2 , we get By substituting the expression of  2 given by (34) into (32) when  = 3 and solving for  3 , we get Next for  ≥ 4, repeat the process of substituting the expression of  (−1) given by into (32) and solving for   to get When  =  is large enough, we may set all the  (+1) in (37) to be zero and obtain Substituting (38) into (37) when  =  − 1, we get Similarly, for  =  − 2,  − 3, . . ., 1, we may perform the substitution of ( (+1) |  0 * * ) into (37) and obtain When  = 1, (40) yields ( 1 |  0 * * ).By using the results given by ( 1 |  0 * * ) and ( 9) to (11), we get the following system of  equations: An inspection of (41) reveals that, among the  equations, only  − 1 of them are linearly independent.Hence, we need to include another linearly independent equation so that the resulting system of  equations has a unique solution.
Equating the sum of the left sides of the equations given by (40) to the sum of the right sides of (40), we get an equation of the form, where the   are constants.
As ∑  =0 ∑  ∑    ≅ 1, we get from (42) an equation involving only  00 , 1 ≤  ≤ .The equation derived from (42) and  − 1 equations chosen from (41) constitute a system of  equations which can be solved to yield numerical answers for  00 , 1 ≤  ≤ .Then using (40), we can get numerical answers for   where  ≥ 1, 0 ≤  ≤ , and 0 ≤  ≤ .The stationary probability that the queue length is  can then be obtained as

Waiting Time Distribution
Suppose a customer arrives at the system which is in the stationary state.Let the time of arrival of the customer be denoted as  = 0. Furthermore, let   be the time the customer needs to wait before being served and the cumulative distribution function (cdf) of the waiting time   .
To find the waiting time distribution   (), we first note that when the system is in the stationary state, an arrival of a customer at time  = 0 which is inside an interval  of length Δ may occur with an approximate probability  +1 Δ if the arrival process is in state  at the beginning of the interval .Meanwhile at the beginning of , the service process may be in state  where 0 ≤  ≤ .Thus the probability that (i) the queue length at the beginning of  is , (ii) the service process is in state  at the beginning of , (iii) the arrival process is in state  at the beginning of , (iv) a customer arrives in  is given approximately by Let  +1 () be the pdf of the service time given that the service process is in state  at the beginning of , and  (−1) () the ( − 1)-fold convolution of ().The customer who arrives in  (see (iv)) under the conditions given by (i), (ii), and (iii) above will have a waiting time of zero if  = 0, and a waiting time in which the pdf is given by the convolution  +1 () *  (−1) () if  ≥ 1.The cdf   () is then given approximately by The cdf   () may also be computed approximately by a simulation procedure described below.Suppose a customer arrives at time  = 0 and the next th customer arrives at time  = ∑  =1   , where  1 ,  2 , . . .are independent and identically distributed with pdf ().Next let the service time of the next th customer be   of which  0 ,  1 , . . .are independent and identically distributed with pdf ().For a chosen large integer , the value of V ̃= {(0, 0 ), ( 1 ,  1 ), . . ., (  ,   )} is generated and the following waiting times are obtained: where  , is the waiting time of the th customer.Then   () ≅ Number of the  , which are less than   . (48)

Numerical Examples
Let Gamma(, ) denote a gamma distribution of which  is the shape parameter and  is the scale parameter.The related probability density function is then given by (; , ) = ( −1  −/ )/(  Γ()).Consider an example in which the service time () has a gamma distribution with the parameter vector ( 1 ,  1 ) = (1.5, 2), and the interarrival time () has another gamma distribution with the parameter vector ( 2 ,  2 ) = (3.1, 2).The utilization factor will then be  = ()/() = 0.48.The reason for considering gamma distribution (, ) with fractional values of the shape parameter  is that the term  −1 appearing in the pdf () and () will usually make the existing analytical methods for finding queue length distribution fail.The reason behind such failure is that when we set  =  + , the function ( + ) −1 cannot be expressed as a finite sum of products of a function of  alone and a function of  alone.
The following is a procedure to find the values of Δ and  (or ).Initially we find the value of  such that the rates at time  ≥  exhibit small variations.A small fractional value (e.g., 0.1 or 0.05) is assigned to Δ and  is then obtained as the integer which is approximately equal to /Δ.If  is very large (e.g.,  > 1000), then a bigger unit is chosen for  until  ≤ 1000.It can be shown that when Δ = 0.04, suitable values of  and  are, respectively,  = 550 and  = 550.By using the proposed numerical method, the stationary queue length distribution is found.The stationary queue length distribution may also be computed using the simulation procedure in the software "QtsPlus" (accompanying software for Gross and Harris [15]) when the number of runs is  1 = 107.The results obtained are shown in Table 1.
From Table 1, we see that the stationary queue length distribution obtained using the proposed numerical method is close to that obtained from the software "QtsPlus." Table 2 shows the stationary queue length probabilities found by the numerical method using  = 30.The reason for choosing  = 30 is that, for each given value of Δ, the value of  30 will then be of the order of about 10 −20 which is small enough for us to get, for 0 ≤  ≤ 7, a value of   which is accurate up to the 6 decimal points.For a given value  of the queue length, we may use the values in the columns under Δ and  in Table 2 to fit a polynomial of low degree in Δ to the queue length probability.The value given by the fitted polynomial when Δ = 0 will represent the final result based on the numerical method for the queue length probability.It can be shown that, for all the values of  considered, the fitted polynomials are all very satisfactory.Thus the final results at Δ = 0 would be quite accurate.Table 2 shows that the numerical results for queue length  = 0, 2, 3, and 4 agree quite well with those based on "QtsPlus." Meanwhile the numerical results for queue length  = 1, 5, 6, and 7 are somewhat more accurate than those based on simulation.
Table 3 shows that the stationary waiting time distribution obtained by using the numerical method in Section 4 is close to that obtained by the simulation procedure.

Discrete Time GI/G/1 Queue
The stationary queue length and stationary waiting time distributions of a discrete time GI/G/1 can also be found by using the proposed numerical method in Sections 2 to 4 after some modifications of the equations for the stationary probabilities given in Section 2. An explanation of why the above modifications are necessary is as follows.
First we note that the values of the   (or   ) for the discrete service time (or arrival time) distribution are such that all the   (or   ) are zero except for the cases when Δ coincides with the service time (or arrival time) which has a nonzero probability of occurrence.Let the values of such  be denoted by  1 ,  2 , . . .,   .The value of a typical    will be such that    Δ is a constant.This means that when Δ is made very small, the value of    will have to be inflated correspondingly.Thus, if the system is not empty at time , the simultaneous occurrence of the events that (A) a customer arrives within the interval (,  + Δ], and (B) a service is completed within the interval (,  + Δ] may not tend to zero when Δ tends to zero.Thus the equations for stationary probabilities given in Section 2 need to be modified by taking into account of the simultaneous occurrence of events (A) and (B).The modified version of the equations in Section 2 is as follows.
When the queue length is  = 0, the values of the   can be found from ( 9)- (11).When the queue length is  = 1, the expressions for   , 1 ≤  ≤ ,  ≤  are the same as those given by ( 16)-( 18), whereas  0 can be computed from the equations as follows: For  ≥ 3, the values of all the   (except  0 ) can be computed using ( 22)-( 23) and ( 50)-(54), whereas those of  0 can be computed using the following equations: We may solve the above equations by using the proposed numerical method in Section 3 to obtain all the values of   and subsequently the stationary queue length distribution.From the values of the stationary probabilities, we can find the stationary waiting time distribution by using the method proposed in Section 4. The cdf   () for the discrete time GI/G/1 queue is now given approximately by Table 4 shows the results of the stationary queue length distribution computed using the proposed numerical method.The table also shows the results given in [11] where the authors found the stationary queue length distribution from the sojourn time distribution using the distributional Little's law.The functions () and () appearing in Tables 4−6 and 8-9 are, respectively, the probability generating functions of the discrete service time and interarrival time.
From Table 4, we can see that the queue length probabilities obtained by using the proposed numerical method are close to those given in [11].When Δ ≤ 1, the values of the   (or   ) are able to capture all the details of the discrete distribution of the service time (or arrival time).Thus the results given in columns 2 and 3 in Table 4 are identical.
Tables 5 and 6 show the stationary queue length distribution in three other examples of discrete queue.
Tables 5 and 6 show that the results obtained by using the proposed numerical method are very close to those given in [11].
From the stationary probabilities, the stationary waiting time distributions can be obtained using (57).The results obtained are shown in Table 7.
For a customer who arrives at time  = 0, his sojourn time is equal to the sum of his waiting time and service time.Thus, from the waiting time and service time distributions of the incoming customer, we can compute his sojourn time distribution.Tables 8 and 9 show the results of the stationary sojourn time distribution computed using the proposed numerical method and those given in [11].
From Tables 8 and 9, we can see that the stationary sojourn time distributions computed by using the proposed numerical method are very close to those given in [11].

Conclusion
Most of the existing methods in the literature find the queue length distribution in the GI/G/1 queue via the waiting time distribution.On the contrary, the present proposed method finds the queue length distribution directly for the CAR/ CAR/1 queue.The accuracy of the numerical results for the distribution can be greatly improved by an extrapolation process.Furthermore the queue length distribution thus found can later be used to find the waiting time distribution and sojourn time distribution.The main drawback of the proposed method is that we may encounter dimensionality problem when  (or ) is very large.For example, when the distribution of the interarrival time tends to a constant only after a long time, the value of  (or ) is very large.To overcome the drawback, we may first obtain a number of large values of Δ such that the values of  and  are not too large, and then find the stationary probabilities for each value of Δ.The stationary probabilities when Δ → 0 may next be obtained by means of extrapolation on the values of Δ.The method proposed in this paper may also be applied to other queuing models which cannot be represented as continuous-time Markov chains.

Table 7 :
Stationary waiting time distributions computed by using the proposed numerical method (Δ = 1.0).