Stationary queue length of a single-server queue with negative arrivals and non-exponential service time distributions

In this paper, a single-server queue with negative customers is considered. The arrival of a negative customer will remove one positive customer that is being served, if any is present. An alternative approach will be introduced to derive a set of equations which will be solved to obtain the stationary queue length distribution. We assume that the service time distribution tends to a constant asymptotic rate when time t goes to infinity. This assumption will allow for finding the stationary queue length of queueing systems with non-exponential service time distributions. Numerical examples for gamma distributed service time with fractional value of shape parameter will be presented in which the steady-state distribution of queue length with such service time distributions may not be easily computed by most of the existing analytical methods.


Introduction
Queues with negative arrivals, known as G-queues, were first introduced by Gelenbe [1, 2] in modelling of neural networks.Since the introduction by Gelenbe, these types of queues have drawn great attention and studied extensively.The application is then widen to other areas such as communication systems, industrial engineering and manufacturing systems [3][4][5].In practical applications, the term of negative customers is used to resemble viruses, inhibitory signals and orders of demand.Typically, the arrival of a negative customer to the queue will remove one positive customer according to some killing disciplines.There are two simplest forms of killing disciplines: removal of a positive customer at the head (RCH) and removal of a positive customer at the end (RCE).Jain and Sigman [6] initiated a disaster killing strategy in which the arrival of a negative customer will remove all positive customers.Boucherie and Boxma [7] then considered the removal of a random number of positive customers upon the arrival of a negative customer.By assuming Poisson arrival of negative customers, Harrison and Pitel [8] derived expressions to find the stationary queue length and sojourn time distributions for an M/M/1 queue.The same authors [9] then extended the model to M/G/1 queue and expressed explicitly the stability conditions.
In this paper, we consider queueing systems with negative arrivals and RCH killing strategy.Numerical method introduced in Koh et al. [10] will be applied to find the stationary queue length distribution.The service time distribution is assumed to have constant asymptotic rate (CAR) when the time t tends to infinity.With this assumption, many distributions can be taken into consideration since in reality, most of the distributions tend to constant when time t goes to infinity.A set of equations for the stationary probabilities will be derived and the stationary queue length distribution can be obtained by solving the equations.Mean queue length computed by the alternative method will be compared to those obtained by the analytical method in [9] and verified with the simulation results.However, we will only find the stationary queue length distribution using the alternative method introduced in this paper since the derivation of probability generating function (pgf) to find the steady-state distribution using Laplace transform have a complex form of expression for an M/G/1 queue.
The rest of the paper is organized as follows.Section 2 gives a description on the queueing model.The set of equations for the stationary probabilities will be derived in Section 3 and the method to solve for the numerical values will be shown in Section 4. In Section 5, expression to find the mean queue length will be derived using the analytical method in [9].Numerical examples will be discussed in Section 6 and the concluding remarks is presented in the Section 7.

Model description
A first come first serve single server queue with negative arrivals is considered.Upon the arrival of a negative customer, one positive customer at the head will be removed if any is present.The negative customer has no effect to the system when the queue is empty.Interarrival time of both the positive and negative customers are assumed to have exponential distribution with rate  and  , respectively.Service time distribution will have a constant asymptotic rate when the time t tends to infinity.Examples of distributions with constant asymptotic rate are exponential, Erlang, Hyperexponential, Gamma and etc.

Derivative of the stationary probabilities
The time axis is first segmented into equal length of interval and denote τk as the interval (( 1) , ]    k t k t for k = 1, 2, 3,….Let f(t) be the probability density function (pdf) of the service time and ) , 1 ( ) where J is large enough such that lim Suppose that a positive customer arrive at time t = 0 and receive service once enter the queue.The probability that the service is completed in the interval τ1 is 1  t .If the service is not completed in the interval τ1, …, τk-1, the probability that the service will be completed in τk is approximately   k t for k = 2, 3, 4, ….where    k J for k ≥ J. Let ,

 
k k and  k be the states of the positive arrival, negative arrival and service processes respectively at the end of τk.We define the state numbers of these processes as follows: ξk = 0, if there is a positive arrival in τk with the rate  , k ≥ 0.
(2) 1, if no positive arrival in τk, k ≥ 0.  k = 0, if  the service for the positive customer ends in τk, for k ≥ 1; or  there is a negative arrival in τk , k ≥ 1; or  the server is idle in τk.Let nk be the number of positive customers at the end of τk.The vector that represents the number of positive customers and the states of the positive arrival, negative arrival and service processes at the end of τk is k nirj p be the probability that at the end of τk, the number of positive customers is n, the states of the positive arrival, negative arrival and service processes is i, r and j, respectively.We assume that the stationary probabilities exist where If the system is not empty at the end of τk-1 and (g) no positive or negative arrival which yields For instance, if Event (e) described above occurs in τk, then the probability that there is one positive customer at the end of the interval τk and When k → ∞, by finding the combinations of 1 k   and the events that could occur in τk which lead to  k , we will obtain the following equations. 3 When n = 2: For n ≥ 2:

Stationary queue length distribution
To solve the set of equations derived in Section 3, the following notations are first introduced: * n P  : Column vector containing all the stationary probabilities of size n.When n = 1, Equations ( 8) to ( 14) can be represented in the matrix form as follows: 25) In general, for n ≥ 3, we can form the following matrix from Equations ( 11) to ( 14) and ( 18) to (21): The same computation in getting (25) will be performed iteratively to obtain the following matrix equation for n ≥ 3: Let n = N be a number that is large enough such that we can set all    N -1, …, 1, by using the iterative procedure, we will obtain: When n = 0, the following equation can be formed from ( 6) and ( 7): Let n = 1 for Equation (30) and by substituting into (31), we get: , the sum of the right hand side (RHS) of (30) for n = 1, 2, …, N and RHS of (32) will yield an equation of the form: where nirj c are constants.
An inspection on the system of equations in (32) reveals that there is one linearly dependent equation.Hence, by substituting one of the equations in (32) with (33), we can then solve for 01 0 r p for 0,1  r .After computing the values of 01 0 r p , we can easily find the numerical values for all the nirj p for n = 1, 2, …, N from Equation (30).The stationary probabilities that there are n positive customers in the queue are then found by

Mean queue length
When all the values of n p are obtained, the mean queue length can be computed using the following function where q L is the number of positive customer in the queue.
On the other side, from Pollaczek-Khinchine formula, the mean queue length is given by MATEC Web of Conferences 189, 02006 (2018) https://doi.org/10.1051/matecconf/201818902006MEAMT 2018 where L, S and  are the queue length random variable, service time random variable and the utilization factor of the queue respectively.
In Harrison and Pitel [9], the modified service time random variable Snew with pdf ) (t b has been found as follows: Substituting (39) into (36), the mean queue length of the M/G/1 queue with negative arrivals and RCH killing strategy can be found by the expression:

Numerical examples
We first present in tables 1 and 2 the mean queue length of the system and compare the numerical results to those obtained by Equation (40).In table 1, service time for the positive customer is assumed to have exponential distribution with rate µ.Then in table 2, Gamma distributed service time is considered with shape parameter κ and scale parameter β.Simulation is also carried out to verify the results.Tables 1 and 2 show that the numerical results obtained by the alternative method is close to those computed from expression (40).The results have been further confirmed by comparison to those obtained from the simulation procedure.

ℒ{ )
In table 3, numerical examples for stationary queue length distributions will be shown in which the service time has a Gamma distribution with fractional values of shape parameter.Steady-state distribution of queue length with such service and/ or interarrival time distributions may not be easily computed by most of the existing analytical methods.
with parameters: From table 3, it can be seen that numerical results computed using the alternative method are close to that obtained by the simulation procedure.

Conclusion
The alternative method discussed in this paper has been successfully applied to find the stationary queue length distribution for queueing system with negative arrivals.The service time distributions are assumed to have constant asymptotic rates when time t tends to infinity.To find the stationary queue length distribution, the alternative method does not involve derivation of pgf which requires Laplace transform on the pdf and cumulative distribution function of the modified service time random variables which may have a complex form when the distribution is general.Simple iterative procedure is used to solve for the stationary probabilities derived in this paper.The drawback of the alternative method is that we may encounter dimensionality problem when J is large.

.
then one of the following events can occur in the next time interval τk:(a) one positive customer enter the queue with the rate  and one negative customer enter the queue with the rate γ and  no positive or negative arrival and no completion of service which yields   that the system is empty at the end of τk-1 and Then only one of the following events can occur in τk:(e) one positive customer enter the queue with the rate  and   1, 0, 0, 1   k ;(f) one negative customer enter the queue with the rate γ and   0, 0,1, 0   k ; (28) into (27) for n = N -1 yields the expression: 37) and its corresponding Laplace transform, ℒ{ ) (t b }(s) is derived as From (38), we can find the utilization factor new  and the second moment of the modified 0) and ℒ ′′{ ) (t b }(0) are the first and second moment of the random variable Snew.

Table 1 .
Comparison of mean queue length when

Table 2 .
Comparison of mean queue length when

Table 3 .
Stationary queue length distribution computed using the alternative numerical method and those obtained from the simulation procedure.