# New Approach for Finding Basic Performance Measures of Single Server Queue

Authors: Siew Khew Koh; Ah Hin Pooi; Yi Fei Tan
Journal: Journal of Probability and Statistics
 (2014)

---

## Abstract

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 functionf(t) and the cumulative distribution function F(t) of the interarrival time are such that the rate f(t)/1-F(t) tends to a constant as t→∞, 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.


---

## Body

## 1. 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.When the interarrival and service times are discrete random variables, the GI/G/1 queue is referred to as a discrete-time queue. Several authors determined the steady-state waiting time distribution in the discrete-time GI/G/1 queue by using the Wiener-Hopf factorization (see [1–5]), the polynomial factorization (see [6, 7]), and the matrix-analytic method (see [8, 9]).Steady-state waiting time distribution in the continuous-time 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 matrix-geometric 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 valueT* such that the probability of the given random variable Y larger than T* is negligible, we may approximate the distribution of Y by a CAR distribution by changing the rate of the distribution of Y at time t for all t>T* to be all equal to the rate at time T*.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Δt. 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 Δt.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 Δt, Haßlinger inspected the queue length L and the time n 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 (L,n) is given by αcLPn where α and c are parameters, and Pn is another parameter which depends on n.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]).


## 2. 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 letg(t) be the probability density function (pdf) of the service time and τk the interval ((k-1)Δt,kΔt] for k=0,1,2,…. Furthermore let
(1)μk=g(kΔt)∫kΔt∞g(u)du,1≤k≤I,
where I is large enough such that
(2)μI≅lim⁡k→∞μk.Suppose a service starts at timet=0. Then, for sufficiently small Δt, the probability that the service will be completed in the interval τ1 is approximately μ1Δt, and, given that the service is not completed in τ1,τ2,…,τk-1, the probability that the service will be completed in τk is approximately μkΔt, k=2,3,4,… where μk=μI for k≥I.For the arrival process, letf(t) be the pdf of the arrival time. Furthermore let
(3)λk=f(kΔt)∫kΔt∞f(u)du,1≤k≤J,
where J is large enough such that
(4)λJ≅lim⁡k→∞λk.Suppose a customer has arrived at timet=0. Then the next customer will arrive in the interval τ1 with an approximate probability λ1Δt, and, given that the next customer does not arrive in the interval τ1,τ2,…,τk-1, the probability that the next customer will arrive in τk will be approximately λkΔt for k=2,3,4,… where λk=λJ for k≥J.Given that a service starts at a time inτ0, we may define the state number ξk of the service process at the end of τk as
(5)ξk={0,if(i)k=0;or(ii)theserviceendsinτk,fork≥1;(ii)or(iii)theserverisidleinτk,min⁡(k,I),iftheservicedoesnotendinτk,k≥1.Next, given that a customer arrives at a time inτ0, we may define the state number ψk of the arrival process at the end of τk as
(6)ψk={0,ifk=0orthenextcustomerarrivesinτk,k≥1,min⁡(k,J),ifthenextcustomerdoesnotarriveinτk,k≥1.Letnk be the queue length at the end of τk and hk=(nk,ξk,ψk). We may refer to hk as the vector of characteristics of the queue at the end of τk.LetPnij(k) be the probability that, at the end of τk, the number of customers in the system is n (including the customer that is being served), the service process is in state i, and the arrival process is in state j, where n≥0, i∈{0,1,2,…,I} and j∈{0,1,2,…,J}. Assume that
(7)Pnij=lim⁡k→∞Pnij(k)
exists. To find the Pnij, we first make the following observations.Suppose at the end ofτk-1 the system is not empty, and the service and arrival processes are in state i-1 and j-1 at the end of τk-1, respectively. Then only one of the following events can occur in the next time interval τk:(a)
a customer enters the system with the arrival rateλj*, and, at the end of τk, the vector of characteristics becomes hk=(n+1,i*,0);(b)
a customer leaves the system with the departure rateμi*, and hk=(n-1,0,j*);(c)
no customers enter or leave the system, andhk=(n,i*,j*).Here,i*=min⁡(i,I) and j*=min⁡(j,J). However if the system is empty at the end of τk-1, and the arrival process is in state j-1, then either one of the following events may occur in τk;(d)
a customer enters the system with arrival rateλj*, and hk=(1,0,0);(e)
no customers enter the system, andhk=(0,0,j*).By setting nk-1=1, i-1=j-1=0 and letting event (b) occurs in τk, we get
(8)P001(k)≅P100(k-1)(μ1Δt)(1-λ1Δt).
When k→∞, (8) yields,
(9)P001≅P100(μ1Δt)(1-λ1Δt).
In general, for a given value of hk, we can likewise find the combinations of hk-1 and the event in τk which lead to hk, and obtain an equation similar to (9). Thus we can obtain the following equations:
(10)P00j≅P00(j-1)(1-λjΔt)+∑i=0j-1P1i(j-1)(μi+1Δt)(1-λjΔt)hhhhhhhhhhhhhfor2≤j≤J-1,(11)P00J≅P00(J-1)(1-λJΔt)+P00J(1-λJΔt)+∑i=0J-1P1i(J-1)(μi+1Δt)(1-λJΔt)+∑i=0J-1P1iJ(μi+1Δt)(1-λJΔt)+P1IJ(μIΔt)(1-λJΔt).
When the queue length is n=1,
(12)Pn00≅∑j=1J-1P00j(λj+1Δt)+P00J(λJΔt),(13)Pn01≅∑i=1I-1P(n+1)i0(μi+1Δt)(1-λ1Δt)hhhh+P(n+1)I0(μIΔt)(1-λ1Δt),(14)Pn0j≅∑i=0,i≠j-1I-1P(n+1)i(j-1)(μi+1Δt)(1-λjΔt)hhhhhh+P(n+1)I(j-1)(μIΔt)(1-λjΔt)hhhhhhhhhhhhhh.hhfor2≤j≤J-1,(15)Pn0J≅∑i=0,i≠J-1I-1P(n+1)i(J-1)(μi+1Δt)(1-λJΔt)hhhhhh+P(n+1)I(J-1)(μIΔt)(1-λJΔt)+∑i=0I-1P(n+1)iJ(μi+1Δt)(1-λJΔt)hhhhhh+P(n+1)IJ(μIΔt)(1-λJΔt),(16)Pnij≅Pn(i-1)(j-1)(1-μiΔt)(1-λjΔt)hhhhhhhfori≤j,1≤i≤I-1,hhhhhhhhhhhhhhhh1≤j≤J-1,(17)PniJ≅Pn(i-1)(J-1)(1-μiΔt)(1-λJΔt)+Pn(i-1)J(1-μiΔt)(1-λJΔt)hhhhhhhhhhhhfor1≤i≤I-1,(18)PnIJ≅Pn(I-1)(J-1)(1-μIΔt)(1-λJΔt)+Pn(I-1)J(1-μIΔt)(1-λJΔt)+PnIJ(1-μIΔt)(1-λJΔt).
When the queue length is n=2, the expressions for Pn0j, 1≤j≤J, are the same as those given by (13), (14), and (15). Other Pnij when the queue length is n=2 can be computed from the equations below:
(19)Pni0≅∑j≥i-1J-1P(n-1)(i-1)j(1-μiΔt)(λj+1Δt)hhhhh+P(n-1)(i-1)J(1-μiΔt)(λJΔt)hhhhhhhhhhhhhhhfor1≤i≤I-1,(20)PnI0≅P(n-1)(I-1)(J-1)(1-μIΔt)(λJΔt)+P(n-1)(I-1)J(1-μIΔt)(λJΔt)+P(n-1)IJ(1-μIΔt)(λJΔt),(21)Pnij≅Pn(i-1)(j-1)(1-μiΔt)(1-λjΔt)hhhhhhhhfori≠j,1≤i≤I-1,hhhhhhhhhhhhhhhh1≤j≤J-1,(22)PniJ≅Pn(i-1)(J-1)(1-μiΔt)(1-λJΔt)+Pn(i-1)J(1-μiΔt)(1-λJΔt)hhhhhhhhhhhhfor1≤i≤I-1,(23)PnIj≅Pn(I-1)(j-1)(1-μIΔt)(1-λjΔt)+PnI(j-1)(1-μIΔt)(1-λjΔt)hhhhhhhhhhhhfor1≤j≤J-1,(24)PnIJ≅Pn(I-1)J(1-μIΔt)(1-λJΔt)+PnI(J-1)(1-μIΔt)(1-λJΔt)+PnIJ(1-μIΔt)(1-λJΔt).When the queue length isn≥3, the values of all the Pnij (except Pni0) can be computed using (21) to (24) and (13) to (15), whereas those of Pni0 can be computed using the following equations:
(25)Pn10≅∑j=1J-1P(n-1)0j(1-μ1Δt)(λj+1Δt)hhhh+P(n-1)0J(1-μ1Δt)(λJΔt),(26)Pni0≅∑j=0,j≠i-1J-1P(n-1)(i-1)j(1-μiΔt)(λj+1Δt)hhhhhh+P(n-1)(i-1)J(1-μiΔt)(λJΔt)hhhhhhhhhhhlhhhhhfor2≤i≤I-1,(27)PnI0≅∑j=0,j≠I-1J-1P(n-1)(I-1)j(1-μIΔt)(λj+1Δt)hhhhh+P(n-1)(I-1)J(1-μIΔt)(λJΔt)+∑j=0J-1P(n-1)Ij(1-μIΔt)(λj+1Δt)hhhhh+P(n-1)IJ(1-μIΔt)(λJΔt).


## 3. Stationary Queue Length Distribution

Before solving (9) to (27) in Section 2, to obtain the stationary queue length distribution, we may first let cij, dij, eij, fj and gij be constants and introduce the following notations:(a)
Pn**={Pnij:0≤i≤I,0≤j≤J};(b)
(Pm**,P(m+1)**,P(m+2)**) denotes the set of equations of the form
(28)∑i=0I∑j=0JcijPmij+∑i=0I∑j=0JdijP(m+1)ij+∑i=0I∑j=0JeijP(m+2)ij≅0;(c)
(Pmij∣P0**,P(m+1)**) denotes the equation of the form
(29)Pmij≅∑j=0JfjP00j+∑i=0I∑j=0JgijP(m+1)ij.With the above notations, (12) to (18) in the case when n=1 can be represented as
(30)(P0**,P1**,P2**),
and (19) to (24) together with (13) to (15) in the case when n=2 may be represented as
(31)(P1**,P2**,P3**).
Furthermore (25) to (27) together with (21) to (24) and (13) to (15) in the case when n≥3 may be represented as
(32)(P(n-1)**,Pn**,P(n+1)**).
It can be shown that, from the set of equations given by (30), we can get
(33)(P1ij∣P0**,P2**)for0≤i≤I,0≤j≤J.
By substituting the expression of P1ij given by (33) into (31), and solving for P2ij, we get
(34)(P2ij∣P0**,P3**)for0≤i≤I,0≤j≤J.
By substituting the expression of P2ij given by (34) into (32) when n=3 and solving for P3ij, we get
(35)(P3ij∣P0**,P4**)for0≤i≤I,0≤j≤J.
Next for n≥4, repeat the process of substituting the expression of P(n-1)ij given by
(36)(P(n-1)ij∣P0**,Pn**),
into (32) and solving for Pnij to get
(37)(Pnij∣P0**,P(n+1)**)for0≤i≤I,0≤j≤J.
When n=N is large enough, we may set all the P(n+1)ij in (37) to be zero and obtain
(38)(PNij∣P0**)for0≤i≤I,0≤j≤J.
Substituting (38) into (37) when n=N-1, we get
(39)(P(N-1)ij∣P0**)for0≤i≤I,0≤j≤J.
Similarly, for n=N-2,N-3,…,1, we may perform the substitution of (P(n+1)ij∣P0**) into (37) and obtain
(40)(Pnij∣P0**)for0≤i≤I,0≤j≤J.
When n=1, (40) yields (P1ij∣P0**). By using the results given by (P1ij∣P0**) and (9) to (11), we get the following system of J equations:
(41)(P00j∣P0**)for0≤j≤J.
An inspection of (41) reveals that, among the J equations, only J-1 of them are linearly independent. Hence, we need to include another linearly independent equation so that the resulting system of J 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,
(42)∑n=1N∑i∑jPnij=∑jkjP00j,
where the kj are constants.As∑n=0N∑i∑jPnij≅1, we get from (42) an equation involving only P00j, 1≤j≤J. The equation derived from (42) and J-1 equations chosen from (41) constitute a system of J equations which can be solved to yield numerical answers for P00j, 1≤j≤J. Then using (40), we can get numerical answers for Pnij where n≥1, 0≤i≤I, and 0≤j≤J. The stationary probability that the queue length is n can then be obtained as
(43)Pn=∑i=0I∑j=0JPnij.


## 4. 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 ast=0. Furthermore, let Wq be the time the customer needs to wait before being served and
(44)Wq(t)=P(Wq≤t)
the cumulative distribution function (cdf) of the waiting time Wq.To find the waiting time distributionWq(t), we first note that when the system is in the stationary state, an arrival of a customer at time t=0 which is inside an interval τ of length Δt may occur with an approximate probability λj+1Δt if the arrival process is in state j at the beginning of the interval τ. Meanwhile at the beginning of τ, the service process may be in state i where 0≤i≤I. Thus the probability that(i)
the queue length at the beginning ofτ is n,(ii)
the service process is in statei at the beginning of τ,(iii)
the arrival process is in statej at the beginning of τ,(iv)
a customer arrives inτis given approximately by
(45)Pnijλj+1Δt.
Let gi+1(t) be the pdf of the service time given that the service process is in state i at the beginning of τ, and g(n-1)(t) the (n-1)-fold convolution of g(t). 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 n=0, and a waiting time in which the pdf is given by the convolution gi+1(t)*g(n-1)(t) if n≥1. The cdf Wq(t) is then given approximately by
(46)Wq(t)≅(∑j=1JP00jλj+1Δt+∑n=1N∑i=0I∑j=0JPnijλj+1Δt∫u=0tgi+1(u)*g(n-1)(u)du)×(∑j=1JP00jλj+1Δt+∑n=1N∑i=0I∑j=0JPnijλj+1Δt)-1.
The cdf Wq(t) may also be computed approximately by a simulation procedure described below.Suppose a customer arrives at timet=0 and the next mth customer arrives at time t=∑k=1mAk, where A1,A2,… are independent and identically distributed with pdf f(t). Next let the service time of the next mth customer be Bm of which B0,B1,… are independent and identically distributed with pdf g(t). For a chosen large integer M, the value of v~={(0,B0),(A1,B1),…,(AM,BM)} is generated and the following waiting times are obtained:
(47)Wq,0=0,Wq,m={0,ifWq,m-1+Bm-1<Am,Wq,m-1+Bm-1-Am,ifWq,m-1+Bm-1≥Am,1≤m≤M-1,
where Wq,m is the waiting time of the mth customer. Then
(48)Wq(t)≅NumberoftheWq,mwhicharelessthantM.


## 5. Numerical Examples

LetGamma(κ,θ) denote a gamma distribution of which κ is the shape parameter and θ is the scale parameter. The related probability density function is then given by f(t;κ,θ)=(tκ-1e-t/θ)/(θκΓ(κ)). Consider an example in which the service time (S) has a gamma distribution with the parameter vector (κ1,θ1)=(1.5,2), and the interarrival time (T) has another gamma distribution with the parameter vector (κ2,θ2)=(3.1,2). The utilization factor will then be ρ=E(S)/E(T)=0.48. The reason for considering gamma distribution (κ,θ) with fractional values of the shape parameter κ is that the term tκ-1 appearing in the pdf f(t) and g(t) will usually make the existing analytical methods for finding queue length distribution fail. The reason behind such failure is that when we set t=x+y, the function (x+y)κ-1 cannot be expressed as a finite sum of products of a function of x alone and a function of y alone.The following is a procedure to find the values ofΔt and I (or J). Initially we find the value of T such that the rates at time t≥T exhibit small variations. A small fractional value (e.g., 0.1 or 0.05) is assigned to Δt and I is then obtained as the integer which is approximately equal to T/Δt. If I is very large (e.g., I>1000), then a bigger unit is chosen for t until I≤1000. It can be shown that when Δt=0.04, suitable values of I and J are, respectively, I=550 and J=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 M1=107. The results obtained are shown in Table 1.Table 1
Comparison of stationary queue length distribution computed from the proposed numerical method, and that obtained from the software “QtsPlus” [(κ1, θ1) = (1.5, 2), (κ2, θ2) = (3.1, 2), Δt = 0.04].
Queue length,n
P  (queue length = n)
Numerical method
Simulation (QtsPlus)
0
0.518854
0.516033
1
0.366244
0.366475
2
0.089790
0.091326
3
0.019748
0.020488
4
0.004227
0.004444
5
8.97
E
-
04
9.48
E
-
04
6
1.90
E
-
04
2.15
E
-
04
7
4.01
E
-
05
5.44
E
-
05
8
8.47
E
-
06
1.33
E
-
05
9
1.79
E
-
06
3.43
E
-
06
10
3.78
E
-
07
1.79
E
-
07
⋮
⋮
⋮
20
6.67
E
-
14
0From Table1, we see that the stationary queue length distribution obtained using the proposed numerical method is close to that obtained from the software “QtsPlus.”Table2 shows the stationary queue length probabilities found by the numerical method using N=30. The reason for choosing N=30 is that, for each given value of Δt, the value of P30 will then be of the order of about 10-20 which is small enough for us to get, for 0≤n≤7, a value of Pn which is accurate up to the 6 decimal points. For a given value n of the queue length, we may use the values in the columns under Δt and n in Table 2 to fit a polynomial of low degree in Δt to the queue length probability. The value given by the fitted polynomial when Δt=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 n considered, the fitted polynomials are all very satisfactory. Thus the final results at Δt=0 would be quite accurate. Table 2 shows that the numerical results for queue length n=0,2,3, and 4 agree quite well with those based on “QtsPlus.” Meanwhile the numerical results for queue length n=1,5,6, and 7 are somewhat more accurate than those based on simulation.Table 2
Comparison of stationary queue length probabilities found by the numerical method using various other values ofΔt, the extrapolated values, and those obtained from the software “QtsPlus” (κ1, θ1) = (1.5, 2), (κ2, θ2) = (3.1, 2).
Δt
Queue length
0
1
2
3
4
5
6
7
0.088000
0.521712
0.366372
0.088018
0.018916
0.003952
0.000818
0.000169
3.47808
E
-
05
0.073333
0.520874
0.366316
0.088550
0.019165
0.004034
0.000841
0.000175
3.63184
E
-
05
0.062857
0.520258
0.366284
0.088935
0.019345
0.004093
0.000858
0.000179
3.74598
E
-
05
0.059459
0.520054
0.366276
0.089060
0.019405
0.004113
0.000864
0.000181
3.78380
E
-
05
0.055000
0.519784
0.366266
0.089226
0.019482
0.004139
0.000871
0.000183
3.83000
E
-
05
0.051163
0.519550
0.366258
0.089369
0.019550
0.004161
0.000878
0.000185
3.87789
E
-
05
0.048889
0.519410
0.366255
0.089454
0.019590
0.004174
0.000882
0.000186
3.90411
E
-
05
0.040000
0.518854
0.366244
0.089790
0.019748
0.004227
0.000897
0.000190
4.00846
E
-
05
0.037931
0.518722
0.366243
0.089868
0.019785
0.004239
0.000900
0.000191
4.03318
E
-
05
0.033846
0.518461
0.366241
0.090024
0.019859
0.004264
0.000908
0.000193
4.08250
E
-
05
0.031429
0.518304
0.366241
0.090117
0.019902
0.004279
0.000912
0.000194
4.11199
E
-
05
0.027848
0.518070
0.366241
0.090254
0.019968
0.004300
0.000918
0.000195
4.15611
E
-
05
0.027500
0.518047
0.366242
0.090268
0.019974
0.004302
0.000919
0.000196
4.16043
E
-
05
Fitted value
0.516162
0.366307
0.091340
0.020486
0.004474
0.000969
0.000210
4.51599
E
-
05
Value based on simulation
0.516033
0.366475
0.091326
0.020488
0.004444
0.000948
0.000215
5.43826
E
-
05Table3 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.Table 3
Comparison of stationary waiting time distribution computed, respectively, by using the proposed numerical method and the simulation procedure (Δt = 0.04).
Time,t
P (Wq≤t)
Numerical method
Simulation
0
0.716771
0.717751
0.04
0.720011
0.719061
0.08
0.723321
0.721774
0.12
0.726602
0.725983
0.16
0.729799
0.728072
0.20
0.732959
0.731658
0.24
0.735939
0.734215
0.28
0.738872
0.738393
0.32
0.741849
0.740451
0.36
0.744914
0.742415
⋮
⋮
⋮
20
0.999523
0.999638


## 6. 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 Sections2 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μk (or λk) for the discrete service time (or arrival time) distribution are such that all the μk (or λk) are zero except for the cases when kΔt coincides with the service time (or arrival time) which has a nonzero probability of occurrence. Let the values of such k be denoted by k1,k2,…,kd. The value of a typical λki will be such that λkiΔt is a constant. This means that when Δt is made very small, the value of λki will have to be inflated correspondingly. Thus, if the system is not empty at time t, the simultaneous occurrence of the events that(A)
a customer arrives within the interval(t,t+Δt], and(B)
a service is completed within the interval(t,t+Δt]may not tend to zero when Δt 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 isn=0, the values of the Pnij can be found from (9)–(11). When the queue length is n=1, the expressions for Pnij, 1≤j≤J, i≤j are the same as those given by (16)–(18), whereas Pn0j can be computed from the equations as follows:
(49)Pn00≅∑j=1J-1P00j(λj+1Δt)+P00J(λJΔt)+∑i=0I-1∑j≥iJ-1Pnij(μi+1Δt)(λj+1Δt)+∑i=0I-1PniJ(μi+1Δt)(λJΔt)+PnIJ(μIΔt)(λJΔt),(50)Pn0j≅∑i=0I-1P(n+1)i(j-1)(μi+1Δt)(1-λjΔt)hhhh+P(n+1)I(j-1)(μIΔt)(1-λjΔt)hhhhhhhhhhhhhhfor1≤j≤J-1,(51)Pn0J≅∑i=0I-1P(n+1)i(J-1)(μi+1Δt)(1-λJΔt)hhhh+P(n+1)I(J-1)(μIΔt)(1-λJΔt)+∑i=0I-1P(n+1)iJ(μi+1Δt)(1-λJΔt)hhhhh+P(n+1)IJ(μIΔt)(1-λJΔt).
When n=2, the values of the Pn0j, 1≤j≤J can be computed using (50)-(51), while (19)-(20) can be used to find the values of Pni0, 1≤i≤I. The values of the PniJ and PnIj can be obtained from (22) and (23), respectively. All the other values of Pnij can be computed using the following equations:
(52)Pn00≅∑i=0I-1∑j=0J-1Pnij(μi+1Δt)(λj+1Δt)+∑i=0I-1PniJ(μi+1Δt)(λJΔt)+∑j=0J-1PnIj(μIΔt)(λj+1Δt)+PnIJ(μIΔt)(λJΔt),(53)Pnij≅Pn(i-1)(j-1)(1-μiΔt)(1-λjΔt)for1≤i≤I-1,1≤j≤J-1,(54)PnIJ≅Pn(I-1)(J-1)(1-μIΔt)(1-λJΔt)+Pn(I-1)J(1-μIΔt)(1-λJΔt)+PnI(J-1)(1-μIΔt)(1-λJΔt)+PnIJ(1-μIΔt)(1-λJΔt).
For n≥3, the values of all the Pnij (except Pni0) can be computed using (22)-(23) and (50)–(54), whereas those of Pni0 can be computed using the following equations:
(55)Pni0≅∑j=0J-1P(n-1)(i-1)j(1-μiΔt)(λj+1Δt)hhhh+P(n-1)(i-1)J(1-μiΔt)(λJΔt)hhhhhhhhhhhhhfor1≤i≤I-1,(56)PnI0≅∑j=0J-1P(n-1)(I-1)j(1-μIΔt)(λj+1Δt)hhhh+P(n-1)(I-1)J(1-μIΔt)(λJΔt)+∑j=0J-1P(n-1)Ij(1-μIΔt)(λj+1Δt)hhhhhh+P(n-1)IJ(1-μIΔt)(λJΔt).We may solve the above equations by using the proposed numerical method in Section3 to obtain all the values of Pnij 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 Wq(t) for the discrete time GI/G/1 queue is now given approximately by
(57)Wq(t)≅(∑j=1JP00jλj+1Δt(1-μi+1Δt)∫u=0tgi+1(u)*g(n-1)(u)du+∑n=1N∑i=0I∑j=0JPnij(λj+1Δt)(1-μi+1Δt)hhhhhhhhhhhh×∫u=0tgi+1(u)*g(n-1)(u)du+∑n=1N∑i=0I∑j=0JPnij(λj+1Δt)(μi+1Δt)∫u=0tg(n-1)(u)du)×(∑j=1JP00jλj+1Δt+∑n=1N∑i=0I∑j=0JPnijλj+1Δt)-1.Table4 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 T(z) and S(z) appearing in Tables 4−6 and 8-9 are, respectively, the probability generating functions of the discrete service time and interarrival time.Table 4
Comparison of stationary queue length distribution computed using the proposed numerical method and that given in Kim and Chaudhry [11].
Queue length,n
Example 1
P (Queue length = n)  [T(z)=z/10+3z2/10+2z3/5+z4/5,S(z)=3z/10+3z2/5+z3/10]
Numerical method(Δt = 1.0)
Numerical method(Δt = 0.1)
Kim and Chaudhry [11]
0
0.333333
0.333333
0.333333
1
0.596799
0.596799
0.596799
2
0.067034
0.067034
0.067034
3
0.002728
0.002728
0.002728
4
0.000101
0.000101
0.000101
5
3.81
E
-
06
3.81
E
-
06
0.000004
6
1.43
E
-
07
1.43
E
-
07
0
7
5.38
E
-
09
5.38
E
-
09
0
8
2.02
E
-
10
2.02
E
-
10
0
9
7.60
E
-
12
7.60
E
-
12
0
10
2.86
E
-
13
2.86
E
-
13
0
⋮
⋮
⋮
⋮From Table4, we can see that the queue length probabilities obtained by using the proposed numerical method are close to those given in [11]. When Δt≤1, the values of the μk (or λk) 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.Tables5 and 6 show the stationary queue length distribution in three other examples of discrete queue.Table 5
Comparison of stationary queue length distribution computed using the proposed numerical method and that given in Kim and Chaudhry [11] (Δt = 1.0).
Queue length,n
Example 2
P (queue length = n) [T(z)=z(z/2+1/2)38,S(z)=(z+z2+⋯+z35)/35]
Numerical method
Kim and Chaudhry [11]
0
0.100000
0.100000
1
0.323533
0.323533
2
0.291648
0.291648
3
0.146160
0.146160
4
0.071163
0.071163
5
0.034640
0.034641
6
0.016862
0.016862
7
0.008208
0.008208
8
0.003995
0.003995
9
0.001945
0.001945
10
9.47
E
-
04
9.47
E
-
04
⋮
⋮
⋮Table 6
Comparison of stationary queue length distributions computed using the proposed numerical method and that given in Kim and Chaudhry [11] (Δt = 1.0).
Queue length,n
Example 3
Example 4
P (Queue length = n)[T(z)=z(z/2+1/2)38,S(z)=(z+z2+z3)/3]
P (Queue length = n)[T(z)=z(z/2+1/2)38,S(z)=(z+z2+⋯+z19)/19]
Numerical method
Kim and Chaudhry [11]
Numerical method
Kim and Chaudhry [11]
0
0.900000
0.900000
0.500000
0.500000
1
0.100000
0.100000
0.494811
0.494811
2
2.49
E
-
12
2.49
E
-
12
0.005189
0.005189
3
2.21
E
-
25
2.21
E
-
25
1.53
E
-
08
1.53
E
-
08
4
1.52
E
-
46
0
2.19
E
-
14
2.19
E
-
14
5
3.54
E
-
60
0
3.13
E
-
20
2.75
E
-
20
6
⋮
⋮
⋮
⋮Tables5 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.Table 7
Stationary waiting time distributions computed by using the proposed numerical method (Δt = 1.0).
Time,t
P (Wq≤t)
Example 1
Example 2
Example 3
Example 4
0
0.847762
0.242019
1
0.955516
1
0.964654
0.259177
1
0.972644
2
0.992963
0.276887
1
0.984262
3
0.998538
0.294817
1
0.991529
4
0.999699
0.313157
1
0.995754
5
0.999938
0.331907
1
0.998002
6
0.999987
0.350904
1
0.999113
7
0.999997
0.370697
1
0.999619
8
0.999999
0.390268
1
0.999842
9
1
0.410278
1
0.999933
⋮
⋮
⋮
⋮
⋮Table 8
Comparison of stationary sojourn time distributions computed by using the proposed numerical method and that given in Kim and Chaudhry [11] (Δt = 1.0).
Time,t
Example 1
Example 2
Sojourn Time Distribution[T(z)=z/10+3z2/10+2z3/5+z4/5,S(z)=3z/10+3z2/5+z3/10]
Sojourn Time Distribution[T(z)=z(z/2+1/2)38,S(z)=(z+z2+⋯+z35)/35]
Numerical method
Kim and Chaudhry [11]
Numerical method
Kim and Chaudhry [11]
0
0
0
0
0
1
0.254329
0.254329
0.006915
0.006915
2
0.543725
0.543708
0.007410
0.007408
3
0.163404
0.163375
0.007911
0.007913
4
0.030347
0.030357
0.008432
0.008429
5
0.006524
0.006540
0.008968
0.008956
6
0.001326
0.001339
0.009503
0.009493
7
2.74
E
-
04
0.000279
0.010059
0.010039
8
5.66
E
-
05
0.000058
0.010610
0.010595
9
1.16
E
-
05
0.000012
0.011178
0.011160
10
2.40
E
-
06
0.000002
0.011746
0.011731
⋮
⋮
⋮
⋮
⋮Table 9
Comparison of stationary sojourn time distributions computed using the proposed numerical method and that given in Kim and Chaudhry [11] (Δt = 1.0).
Time,t
Example 3
Example 4
Sojourn Time Distribution[T(z)=z(z/2+1/2)38,S(z)=(z+z2+z3)/3]
Sojourn Time Distribution[T(z)=z(z/2+1/2)38,S(z)=(z+z2+⋯+z19)/19]
Numerical method
Kim and Chaudhry [11]
Numerical method
Kim and Chaudhry [11]
0
0
0
0
0
1
0.333333
0.333333
0.050290
0.050290
2
0.333333
0.333333
0.051192
0.051193
3
0.333333
0.333333
0.051803
0.051804
4
1.62
E
-
11
1.62
E
-
11
0.052186
0.052187
5
3.97
E
-
13
4.04
E
-
13
0.052408
0.052408
6
⋮
⋮
⋮
⋮For a customer who arrives at timet=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 Tables8 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].


## 7. 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 whenI (or J) is very large. For example, when the distribution of the interarrival time tends to a constant only after a long time, the value of I (or J) is very large. To overcome the drawback, we may first obtain a number of large values of Δt such that the values of I and J are not too large, and then find the stationary probabilities for each value of Δt. The stationary probabilities when Δt→0 may next be obtained by means of extrapolation on the values of Δt.The method proposed in this paper may also be applied to other queuing models which cannot be represented as continuous-time Markov chains.



---

*Source: 851738-2014-10-20.xml*