The Geo/G^a,Y/1/N Queue Revisited

Abstract

We not only present an alternative and simpler approach to find steady-state distributions of the number of jobs for the finite-space queueing model $Geo/G^{a,Y}/1/N$ using roots of the inherent characteristic equation, but also correct errors in some published papers. The server has a random service capacity Y, and it processes the jobs only when the number of jobs in the system is at least ‘a’, a threshold value. The main advantage of this alternative process is that it gives a unified approach in dealing with both finite- and infinite-buffer systems. The queue-length distribution is obtained both at departure and random epochs. We derive the relation between the discrete-time Geo/$G^{a,Y}$/1/N queue and its continuous-time analogue. Finally, we deal with performance measures and numerical results.

1. Introduction

Discrete-time queues with batch service have numerous applications in various areas such as transportation systems, traffic, production, manufacturing, telecommunication, and cloud computing. In various real-life settings, it is often noted that the jobs are served in batches. The server may serve with fixed maximum or variable capacity in batch service systems. For more details on batch service queues, one may refer to Chaudhry and Templeton [1] as well as Medhi [2]. Discrete-time queues are more notable in systems’ modelling, see [3,4,5].

In discrete-time queues, it is assumed that arrivals and departures occur at boundary epochs of time slots. Further, discrete-time queues deal with an early arrival system (EAS) or a late arrival system with delayed access (LAS-DA). For more on this, see Hunter [4]. We may note that EAS and LAS-DA policies are similar to departure-first (DF) and arrival-first (AF), respectively; see Gravey and Hébuterne [6].

Several researchers study single-server batch-service discrete-time queues with various phenomena, see [7,8,9,10,11,12,13]. In Gupta and Goswami [10], they discuss a discrete-time finite-buffer general bulk service queue under both LAS-DA and EAS policies. The model involving batch-size-dependent service in a discrete-time queue where inter-arrival times and the service times follow geometric and general distribution, respectively, has been discussed by Banerjee et al. [14]. In Yi et al. [13], the authors analyze a discrete-time Geo/$G^{a,Y}$/1/N queue, where service is initiated only when the number of jobs in the system is at least ‘a’. In Zeng and Xia [15], the authors discuss M/$G^{a,b}$/1/N queue where service is in batches with minimum threshold a, maximum capacity b and the buffer size, N, finite or infinite.

At some point, finding the roots of the characteristic equation seemed difficult, mainly when the number involved was large. Several researchers have made these comments. Given this, the procedure for solving queueing models led to the matrix-analytic or matrix-geometric method. In this connection, see the remark below. Following Chaudhry et al. [16], other researchers have attempted to show that roots can be found, and using them leads to nice analytic and numerical solutions. Gouweleeuw [17] shows that the roots’ approach for finding probabilities from generating functions in precise expressions is effective. Further, in the case of large buffer size, solving simultaneous equations gives rise to poor reliability and takes considerable time; see Powell [18] (p. 141), who, while dealing with the model $M/D^c/1$, states that when $c\ge 40$ the method using simultaneous equations breaks down leading to negative probabilities. The goal of this paper is to give an alternative solution that is analytically powerful, simple, and easy to implement numerically. It may be stated further that this method has not been used to discuss the discrete-time queueing system that deals with batch services and a finite buffer.

In real-world systems, we encounter many finite-buffer systems such as telecommunication networks. Because of this, we study the $Geo/G^{a,Y}/1/N$ model under the assumption of late arrival and delayed access system (LAS-DA). Here, we assume that the single server with variable service capacity will process the jobs only when there are at least $‘a^{’}$ jobs in the system. In Yi et al. [13], the authors found the queue-length distribution at post-departure by solving simultaneous equations and random epoch by applying the “rate in = rate out” arguments. We develop an alternate process to find the queue-length distributions at post-departure and random epochs.

The principal contributions of this work may be summed up as follows:

• We find an alternative method to obtain the steady-state queue-length distributions of $Geo/G^{a,Y}/1/N$ at post-departure and random epochs.
• The approach presented in this paper unifies in a way that can handle both the infinite-space as well as finite-space models at the same time.
• We point out the incorrectness of queue-length distributions’ numerical results (at random epochs) reported in Yi et al. [13]. They also assumed batches with a random capacity Y having probability mass function (pmf) $P(Y=i)=y_i,i=0,1,\dots ,b$ instead of $i=a,a+1,\dots ,b$.
• We compute the steady-state queue-length distributions of $Geo/G^{a,Y}/1/N$ at post-departure and random epochs when $\rho >1$, which is missing in Yi et al. [13]. Further, we point out the incorrectness of the formula for the mean waiting time in the queue (using Little’s rule) in Yi et al. [13].
• We can obtain the continuous-time solution for the model $M/G^{a,Y}/1/N$ (see Appendix A) and the procedure used here can be applied to obtain a solution for this continuous-time model too. Further, it is anticipated that, using this method, we can obtain waiting-time distribution using Little’s law, a problem for which no solution is available, even using the matrix-analytic method. The primary purpose of this paper is to show its unifying power and superiority over other methods, and to give a simple solution to the existing problem.
• Finally, we compare the roots’ method against the process that uses simultaneous equations and present the results in the numerical section. It clearly shows that the roots approach takes less time.

The remaining paper is structured as follows. Section 2 specifies the model. Section 3 analyzes the $Geo/G^{a,Y}/1/N$ system and finds queue-length distributions for the LAS-DA policy. Section 4 examines various system performance measurements. Section 5 provides some numerical results and, finally, the paper is concluded in Section 6.

Remark 1.

It may be useful to comment on the matrix-analytic and the root-finding method. Kendall [19] made a famous remark that queueing theory wears the Laplacian curtain. Kleinrock [20] (p. 291) states, “One of the most difficult parts of this method of spectrum factorization is to solve for the roots”. Neuts (see Neuts’ book [21] and also Stidham [22]) states, “In discussing matrix-analytic solutions, I had pointed out that when the Rouché roots coincide or are close together, the method of roots could be numerically inaccurate. When I finally got copies of Crommelin’s papers, I was elated to read that, for the same reasons as I, he was concerned about the clustering of roots. In 1932, Crommelin knew; in 1980, many of my colleagues did not …”. Following this, several other researchers made similar comments. Given this, the procedure for solving queueing models led to the matrix-analytic or matrix-geometric method. In response, Chaudhry et al. [16] showed that the roots can be found even when the number involved is large. (This was done when MATHEMATICA OR MAPLE failed to give those roots. We can now use this software to find roots.)

2. Model Description

We consider a $Geo/G^{a,Y}/1/N$ queue where jobs arrive following a Bernoulli process with parameter $\lambda$. A single server processes the jobs on a first-come-first-served (FCFS) discipline in batches with a random capacity Y possessing a probability mass function (pmf) $P(Y=i)=y_i,i=a,a+1,\dots ,b$ with ${\sum }_{i=a}^b{y}_i=1$, the probability generating function (PGF) $Y\left(z\right)={\sum }_{i=a}^b{y}_i{z}^i$, mean $E\left(Y\right)={\sum }_{i=a}^{b}i{y}_i$ and $Y^{{}^{\prime }}\left(z\right)=z^{b}Y\left(z^{-1}\right)$. We assume the minimum and maximum threshold values of the random variable Y as a and b, respectively. When there are at least $‘a{}^{’}$ jobs in the queue, the server commences serving a batch of size i with probability $y_i$ (when there are $Y\le b$, it takes all of them). When the number of jobs comes down below a threshold value $a(\ge 1)$ in the system, the server remains idle but awaits the number of jobs to rise to a; when it attains a, it resumes service. The service times $\{S_n,n\ge 1\}$ are independently and identically distributed (iid) with arbitrary pmf $s_k=P(S_n=k)$, $k=1,2,\dots$ and $s_0=0$, the PGF $S\left(z\right)={\sum }_{i=1}^{\infty }{s}_i{z}^i$ and mean service time $E\left(S\right)=s=\frac{d}{dz}{S\left(z\right)|}_{z=1}=1/\mu$.

The processing times of the jobs are independent of the arrival process and the number of jobs served. The waiting buffer has a finite capacity N with $b\le N$. Thus, in the system, no more than $N+b$ jobs can be available anytime. We presume offered load of the system as $\rho =\lambda E\left(S\right)/E\left(Y\right)$. In LAS-DA policy, arrivals occur in $(u-,u)$, and departures take place in $(u,u+)$; arrivals supersede departures. Figure 1 describes the various time periods at which events occur. For more details on this, one may refer to [5,6].

3. Queue-Length Distributions

Here, we find steady-state queue-length distributions at various epochs of $Geo/G^{a,Y}/1/N$ queue with the LAS-DA policy.

3.1. Post-Departure Epoch Probabilities

Let $Q_{\ell }^{+}$ be the number of jobs in the queue after completing the ℓth service. Suppose $A_{\ell +1}$ and $Y_{\ell +1}$ represent the number of arrivals throughout the processing time on the $(\ell +1)$th job and the processing capacity of $(\ell +1)$th service, respectively. As per the batch service rule, the departure epoch queue lengthis

$$Q_{\ell +1}^{+}=\min \left({\left(Q_{\ell }^{+}-Y_{\ell +1}\right)}^{+}+A_{\ell +1},N\right),$$

where $x^{+}$ = max $(x,0)$. The probability distribution of $A_{\ell +1}$ is

$$\begin{array}{cc}\hfill k_n& =P(A_{\ell +1}=n)={\displaystyle \sum _{j=1}^{\infty }}P\left(A_{\ell +1}=n|S_{\ell +1}=j\right)P\left(S_{\ell +1}=j\right)\hfill \\ \hfill & ={\displaystyle \sum _{j=n}^{\infty }}{s}_j\left(\genfrac{}{}{0pt}{}{j}{n}\right){\lambda }^n{(1-\lambda )}^{j-n},n\ge 0.\hfill \end{array}$$

Here we may note that arrivals are generated by a Bernoulli sequence by the property of geometric interarrival times. In LAS-DA, if the service time of the $(\ell +1)$th job is j slots, then there will be j time slots where n arrivals may occur. One may note that $\left(\genfrac{}{}{0pt}{}{j}{n}\right){\lambda }^n{(1-\lambda )}^{j-n}$ is the probability of n arrivals in j slots. Let $K\left(z\right)={\sum }_{n=0}^{\infty }{k}_n{z}^n$ be the probability generating function of the sequence $\{k_n,n=0,1,\dots \}$. Thus

$$\begin{array}{cc}\hfill K\left(z\right)& ={\displaystyle \sum _{n=0}^{\infty }}{\displaystyle \sum _{j=n}^{\infty }}{s}_j\left(\genfrac{}{}{0pt}{}{j}{n}\right){\lambda }^n{(1-\lambda )}^{j-n}{z}^n={\displaystyle \sum _{j=0}^{\infty }}{s}_j{(1-\lambda +\lambda z)}^j=S(1-\lambda +\lambda z).\hfill \end{array}$$

Transition probabilities in one step of underlying Markov chain for $p_{ij}=\underset{\ell \to \infty }{lim}Pr\{Q_{\ell +1}^{+}=j|Q_{\ell }^{+}=i\}$ are given as

$$p_{ij}=\left\{\begin{array}{cc}{k}_j,\hfill & 0\le i\le a,0\le j\le N-1,\hfill \\ k_j{\displaystyle \sum _{r=i}^b}{y}_r+{\displaystyle \sum _{r=a}^{i-1}}{y}_r{k}_{j-i+r},\hfill & a+1\le i\le b-1,i-a-1\le j\le N-1,\hfill \\ {\displaystyle \sum _{r=a}^b}{y}_r{k}_{j-i+r},\hfill & b\le i\le N,i-b\le j\le N-1,\hfill \\ {\ell }_N,\hfill & 0\le i\le a,j=N,\hfill \\ {\ell }_N{\displaystyle \sum _{r=i}^b}{y}_r+{\displaystyle \sum _{r=a}^{i-1}}{y}_r{\ell }_{N-i+r},\hfill & a+1\le i\le b-1,j=N,\hfill \\ {\displaystyle \sum _{r=a}^b}{y}_r{\ell }_{N-i+r},\hfill & b\le i\le N,j=N,\hfill \end{array}\right.$$

(1)

where ${\ell }_j={\displaystyle \sum _{r=j}^{\infty }}{k}_r$ and $k_j$ with $j<0$ defined to be zero, and which leads to the transition probability matrix $P=\left(p_{ij}\right)$ as

$$P=\begin{array}{cc}& \begin{array}{ccccccc}0& 1& \dots & j& \dots & N-1& N\end{array}\\ \begin{array}{c}0\\ 1\\ ⋮\\ a\\ a+1\\ & ⋮\\ b\\ b+1\\ ⋮\\ N\end{array}& \left(\begin{array}{ccccccc}{k}_0& k_1& \dots & k_j& \dots & k_{N-1}& {\ell }_N\\ k_0& k_1& \dots & k_j& \dots & k_{N-1}& {\ell }_N\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ k_0& k_1& ⋮& k_j& \dots & k_{N-1}& {\ell }_N\\ k_0{\displaystyle \sum _{r=a+1}^b}{y}_r& k_1{\displaystyle \sum _{r=a+1}^b}{y}_r& ⋮& k_j{\displaystyle \sum _{r=a+1}^b}{y}_r& \dots & k_{N-1}{\displaystyle \sum _{r=a+1}^b}{y}_r& {\ell }_N{\displaystyle \sum _{r=a+1}^b}{y}_r\\ & +k_0{y}_a& & +k_{j-1}{y}_a& & +k_{N-2}{y}_a& +{\ell }_{N-1}{y}_a\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ k_0{y}_b& {\displaystyle \sum _{r=b-1}^b}{y}_r{k}_{1-b+r}& ⋮& {\displaystyle \sum _{r=a}^b}{y}_r{k}_{j-b+r}& \dots & {\displaystyle \sum _{r=a}^b}{y}_r{k}_{N-1-b+r}& {\displaystyle \sum _{r=a}^b}{y}_r{k}_{N-b+r}\\ 0& k_0{y}_b& ⋮& {\displaystyle \sum _{r=a}^b}{y}_r{k}_{j-b+r-1}& \dots & {\displaystyle \sum _{r=a}^b}{y}_r{k}_{N-b+r-1}& {\displaystyle \sum _{r=a}^b}{y}_r{k}_{N-b+r-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& {\displaystyle \sum _{r=a}^b}{y}_r{k}_{j-N+r}& ⋮& {\displaystyle \sum _{r=a}^b}{y}_r{k}_{r-1}& {\displaystyle \sum _{r=a}^b}{y}_r{\ell }_r\end{array}\right)\end{array}.$$

The $(\ell +1)$th batch service starts when there are $‘a{}^{’}$ jobs in the queue, and the state changeover occurs from $0\le i\le a-1$ to $0\le j\le N-1$. If the state changeover is from state $i\ge a$ to state $0\le j\le N-1$, then there is a busy period between the leaving epoch of the ℓth batch and the commencement of processing $(\ell +1)$st batch. Let the steady-state probability $p^{+}=\{P_0^{+},P_1^{+},\dots ,P_n^{+},\dots ,P_N^{+}\}$ represent ℓ jobs at departure epochs. Then, in steady-state, $p^{+}=p^{+}P$ can be expressed as follows:

$$\begin{array}{cc}\hfill P_j^{+}& ={\displaystyle \sum _{i=0}^a}{P}_i^{+}{k}_j+{\displaystyle \sum _{i=a+1}^{b-1}}{P}_i^{+}{\displaystyle \sum _{r=i}^b}{y}_r{k}_j+{\displaystyle \sum _{i=a+1}^{b-1}}{P}_i^{+}{\displaystyle \sum _{r=a}^{i-1}}{y}_r{k}_{j-i+r}+{\displaystyle \sum _{i=b}^N}{P}_i^{+}{\displaystyle \sum _{r=a}^b}{y}_r{k}_{j-i+r},\hfill \\ & 0\le j\le N-1,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill P_N^{+}& ={\displaystyle \sum _{i=0}^a}{P}_i^{+}{\ell }_N+{\displaystyle \sum _{i=a+1}^{b-1}}{P}_i^{+}{\displaystyle \sum _{r=i}^b}{y}_r{\ell }_N+{\displaystyle \sum _{i=a+1}^{b-1}}{P}_i^{+}{\displaystyle \sum _{r=a}^{i-1}}{y}_r{\ell }_{N-i+r}+{\displaystyle \sum _{i=b}^N}{P}_i^{+}{\displaystyle \sum _{r=a}^b}{y}_r{\ell }_{N-i+r},\hfill \end{array}$$

(3)

where the normalization condition is ${\displaystyle \sum _{j=0}^N}{P}_j^{+}=1$. It may be noted that Equation (3) is redundant and will not be considered in analysis hereafter. Specify PGF of $P_j^{+}$ as $P^{+}\left(z\right)={\displaystyle \sum _{j=0}^N}{P}_j^{+}{z}^j$. Multiplying Equation (2) by $z^j$ and then adding overall j, we obtain

$$\begin{array}{cc}\hfill P^{+}\left(z\right)& ={\displaystyle \sum _{j=0}^{N-1}}{k}_j{z}^j{\displaystyle \sum _{i=0}^a}{P}_i^{+}+{\displaystyle \sum _{j=0}^{N-1}}{z}^j{\displaystyle \sum _{i=a+1}^{b-1}}{P}_i^{+}\left({\displaystyle \sum _{r=i}^b}{y}_r\right)k_j\hfill \\ & +{\displaystyle \sum _{j=0}^{N-1}}{z}^j{\displaystyle \sum _{i=a+1}^{b-1}}{P}_i^{+}{\displaystyle \sum _{r=a}^{i-1}}{y}_r{k}_{j-i+r}+{\displaystyle \sum _{j=0}^{N-1}}{z}^j{\displaystyle \sum _{i=b}^N}{P}_i^{+}{\displaystyle \sum _{r=a}^b}{y}_r{k}_{j-i+r}+P_N^{+}{z}^N,\hfill \end{array}$$

$$\begin{array}{c}{P}^{+}\left(z\right)\left[1-K\left(z\right)Y\left(\frac{1}{z}\right)\right]=K\left(z\right){\displaystyle \sum _{i=0}^{a-1}}{P}_i^{+}\left(1-z^{i}Y\left(\frac{1}{z}\right)\right)+P_N^{+}{z}^N\hfill \\ +K\left(z\right){\displaystyle \sum _{i=a}^b}{P}_i^{+}\left({\displaystyle \sum _{r=i}^b}{y}_r-z^i{\displaystyle \sum _{r=i}^b}\frac{y_r}{z^r}\right)-{\displaystyle \sum _{i=a}^b}{y}_i{\displaystyle \sum _{j=i+1}^N}{P}_j^{+}{\displaystyle \sum _{r=N-j+a}^{\infty }}{k}_r{z}^{r+j-i}\hfill \\ -\left({\displaystyle \sum _{i=0}^{a-1}}{P}_i^{+}+{\displaystyle \sum _{i=a}^b}{P}_i^{+}{\displaystyle \sum _{r=i}^b}{y}_r\right){\displaystyle \sum _{j=N}^{\infty }}{k}_j{z}^j.\hfill \end{array}$$

Simplifying the above equation, we get

$$\begin{array}{cc}\hfill P^{+}\left(z\right)& =\frac{K\left(z\right)\left[{\displaystyle \sum _{i=0}^{a-1}}{P}_i^{+}\left(z^b-z^i{Y}^{{}^{\prime }}\left(z\right)\right)+{\displaystyle \sum _{i=a}^{b-1}}{P}_i^{+}\left({\displaystyle \sum _{r=i}^b}{y}_r{z}^b-z^i{\displaystyle \sum _{r=i}^b}{y}_r{z}^{b-r}\right)\right]}{z^b-K\left(z\right)Y^{{}^{\prime }}\left(z\right)}\hfill \\ \hfill & +\frac{z^{N+b}\left(P_N^{+}-{\displaystyle \sum _{i=a}^b}{y}_i{\displaystyle \sum _{j=i+1}^N}{P}_j^{+}{\displaystyle \sum _{r=N-j+a}^{\infty }}{k}_r{z}^{r+j-i-N}\right)}{z^b-K\left(z\right)Y^{{}^{\prime }}\left(z\right)}\hfill \\ \hfill & -\frac{z^{N+b}\left({\displaystyle \sum _{i=0}^{a-1}}{P}_i^{+}+{\displaystyle \sum _{i=a}^b}{P}_i^{+}{\displaystyle \sum _{r=i}^b}{y}_r\right){\displaystyle \sum _{j=N}^{\infty }}{k}_j{z}^{j-N}}{z^b-K\left(z\right)Y^{{}^{\prime }}\left(z\right)}.\hfill \end{array}$$

(4)

Only the first expression on the right side of the Equation (4) will add to the coefficient of $z^j,j=0,1,\dots N$. To the right of Equation (4), we ignore the second and third expressions, which consist of an output of z higher than $N+b$. These are not required as we want to compare the coefficients of $z^j$ for $j\le N$ on both sides in Equation (4) to find $P_j^{+}$ for $j=0,1,\dots ,N$. Let

$$\begin{array}{c}\hfill P_N^{+}\left(z\right)=\frac{K\left(z\right)\left[{\displaystyle \sum _{i=0}^{a-1}}{P}_i^{+}\left(z^b-z^i{Y}^{{}^{\prime }}\left(z\right)\right)+{\displaystyle \sum _{i=a}^{b-1}}{P}_i^{+}\left({\displaystyle \sum _{r=i}^b}{y}_r{z}^b-z^i{\displaystyle \sum _{r=i}^b}{y}_r{z}^{b-r}\right)\right]}{z^b-K\left(z\right)Y^{{}^{\prime }}\left(z\right)},\end{array}$$

(5)

which is equivalent to the PGF of an infinite buffer case. The function $P_N^{+}\left(z\right)$ is fully determined once $P_i^{+},i=0,1,\dots ,b-1$ are known. One may observe that when $\rho <1$, Equation (5) denotes the PGF of discrete-time infinite buffer $Geo/G^{a,Y}/1$ queue. We can calculate the probabilities for $\rho \ge 1$ in the case of a finite buffer $Geo/G^{a,Y}/1/N$ queue.

Remark 2.

Using $a=1$, the model is reduced to $Geo/G^Y/1/N$ queue.

Remark 3.

Taking $y_1=1$, $y_i=0,\forall 2\le i\le b$ and $a=1$, the model becomes $Geo/G/1/N$ queue and Equation (5) establishes PGF as $P_N^{+}\left(z\right)=\frac{P_0^{+}K\left(z\right)\left(1-z\right)}{K\left(z\right)-z}$, where $K\left(z\right)=S(1-\lambda +\lambda z)$, which corresponds to the results of [23].

Remark 4.

Taking $y_b=1$ and $y_i=0,\forall i\ne b$, the model is reduced to $Geo/G^{(a,b)}/1/N$ and Equation (5) establishes the PGF as $P_N^{+}\left(z\right)=\frac{K\left(z\right){\displaystyle \sum _{i=0}^{b-1}}{P}_i^{+}\left(z^b-z^i\right)}{z^b-K\left(z\right)}$.

Intending to establish a unified method to solve the queueing system $Geo/G^{a,Y}/1/N$, we obtain ${\left\{P_n^{+}\right\}}_0^N$ from $P^{+}\left(z\right)$ by using the roots of characteristic equation and partial-fraction expansion. The literature on queueing systems shows that arrival/service-time distributions that possess the generating function as a rational function deal with the broad class of distributions see [24]. For this, we suppose that $K\left(z\right)=S(1-\lambda +\lambda z)$ as a rational function in z, specified by

$$\begin{array}{cc}\hfill K\left(z\right)& =S(1-\lambda +\lambda z)=\frac{f\left(z\right)}{g\left(z\right)},\hfill \end{array}$$

where $f\left(z\right)$ and $g\left(z\right)$ are polynomials of degree m and s, respectively, where m and s can have any value, e.g., m can be greater than s, e.g., see Example 3(ii). Thus, we have from Equation (5),

$$\begin{array}{cc}\hfill P_N^{+}\left(z\right)& =\frac{f\left(z\right)\left[{\displaystyle \sum _{i=0}^{a-1}}{P}_i^{+}\left(z^b-z^i{Y}^{{}^{\prime }}\left(z\right)\right)+{\displaystyle \sum _{i=a}^{b-1}}{P}_i^{+}\left({\displaystyle \sum _{r=i}^b}{y}_r{z}^b-z^i{\displaystyle \sum _{r=i}^b}{y}_r{z}^{b-r}\right)\right]}{z^{b}g\left(z\right)-f\left(z\right)Y^{{}^{\prime }}\left(z\right)}.\hfill \end{array}$$

(6)

The denominator of Equation (6) is a polynomial of degree $b+s$ which when equated to zero has $b+s$ roots inside, on, or outside the unit circle $\mid z\mid =1$, say, ${\gamma }_1=1,{\gamma }_i(i=2,3,\dots ,b+s)$.

Remark 5.

If the denominator $z^{b}g\left(z\right)-f\left(z\right)Y^{{}^{\prime }}\left(z\right)$ of Equation (6) $=0$ obtains roots close to each other or repeated roots, we may obtain them by applying advanced software packages, for instance, MATHEMATICA or MAPLE. The MAPLE script for calculating repeated roots is exemplified below for Equation

$$u\left(y\right)=(y-2){(y-5)}^2{(y-7)}^3(y-11).$$

$$\begin{array}{ccc}& \begin{array}{cc}& restart:Digits:=10:with\left(RootFinding\right):\hfill \\ & m:=\left(y-2\right){\left(y-5\right)}^2{\left(y-7\right)}^3\left(y-11\right);\hfill \\ & Analytic(u,y,re=-1..10,im=-2..10);\hfill \\ & 7.00000000000000,7.00000000000000,7.00000000000000,5.00000000000000,\hfill \\ & 5.00000000000000,2.00000000000000,11.00000000000000\hfill \end{array}\hfill & \end{array}$$

According to Rouché’s theorem, the denominator of Equation (6) has b zeros say, ${\gamma }_i(i=1,2,\dots ,b)$ inside the unit circle. As $P_N^{+}\left(z\right)$ converges in $\mid z\mid \le 1$, the b zeros within the unit circle of the denominator should cancel with the b zeros of the numerator. After canceling the b factors in the numerator and denominator, we can re-write Equation (6) as

$$\begin{array}{cc}\hfill P_N^{+}\left(z\right)& =T\left(C\left(z\right)+\frac{f_1\left(z\right)}{{\displaystyle \prod _{i=b+1}^{b+s}}(z-{\gamma }_i)}\right),\hfill \end{array}$$

(7)

where $C\left(z\right)={\displaystyle \sum _{i=0}^{n_0}}{c}_i{z}^i$ and T is a normalizing constant. Note that when $2m<s$, then $C\left(z\right)$ will be zero. In the partial-fraction process, a slight modification is needed ([25], p. 221) when all the roots are not distinct. Since we are looking at the finite buffer queue system, three instances appear here.

• If $\rho <1$ the s roots ${\gamma }_i,i=b+1,b+2,\dots ,b+s$ remain outside the circle $\mid z\mid =1$.
• If $\rho =1$, among the s roots, one root is $‘1^{’}$, and the other roots ${\gamma }_i,i=b+2,b+3,\dots ,b+s$ are outside the unit circle $\mid z\mid =1$.
• If $\rho >1$, among s roots, one root is inside, say ${\gamma }_{b+1}$ and the other roots ${\gamma }_i,i=b+2,b+3,\dots ,b+s$ are outside, see [26]. One may note that when $\rho >0$ increases, one positive real root comes closer to the origin from right to left.

The expression (7) is tractable for inversion. Applying partial-fraction expansion to Equation (7) yields

$$\begin{array}{cc}\hfill P_N^{+}\left(z\right)& =T\left(C\left(z\right)+{\displaystyle \sum _{i=b+1}^{b+s}}\frac{M_i}{z-{\gamma }_i}\right),\hfill \end{array}$$

(8)

where

$$\begin{array}{cc}\hfill M_i& =\frac{f_1\left({\gamma }_i\right)}{{\displaystyle \prod _{j=b+1,j\ne i}^{b+s}}({\gamma }_i-{\gamma }_j)}.\hfill \end{array}$$

Using Equation (8), we have

$$\begin{array}{cc}\hfill P_n^{+}& =\left\{\begin{array}{cc}T\left(c_n+{\displaystyle \sum _{i=b+1}^{b+s}}\frac{-M_i}{{\gamma }_i^{n+1}}\right),\hfill & \mathrm{if}n=0,1,\cdots ,n_0,\hfill \\ T{\displaystyle \sum _{i=b+1}^{b+s}}\frac{-M_i}{{\gamma }_i^{n+1}},\hfill & \mathrm{if}n=n_0+1,n_0+2,\cdots ,N.\hfill \end{array}\right.\hfill \end{array}$$

(9)

Employing the normalization condition ${\displaystyle \sum _{n=0}^N}{P}_n^{+}=1$, we obtain the only unknown T as

$$\begin{array}{cc}\hfill T& =\left\{\begin{array}{cc}{\left({\displaystyle \sum _{i=0}^{n_0}}{c}_i-{\displaystyle \sum _{i=b+1}^{b+s}}\frac{M_i}{{\gamma }_i}\times \frac{1-{\gamma }_i^{-(N+1)}}{1-{\gamma }_i^{-1}}\right)}^{-1},\hfill & \mathrm{if}\rho \ne 1,\hfill \\ {\left({\displaystyle \sum _{i=0}^{n_0}}{c}_i-M_{b+1}(N+1)-{\displaystyle \sum _{i=b+2}^{b+s}}\frac{M_i}{{\psi }_i}\times \frac{1-{\psi }_i^{-(N+1)}}{1-{\psi }_i^{-1}}\right)}^{-1},\hfill & \mathrm{if}\rho =1.\hfill \end{array}\right.\hfill \end{array}$$

(10)

Thus, once all the roots are known, we can get the distribution for the number in queue.

Remark 6.

It may be noted that it is also possible to find the probabilities ${\left\{P_n^{+}\right\}}_0^N$ by assuming the solution of the form $P_j^{+}=C{z}^j$, where $C\ne 0$. Unfortunately, if we use this method, we have to solve for N simultaneous equations.

3.2. Relationship between the Queue-Length Distributions at Post-Departure and Random Epochs

This sub-section establishes associations between probability at random and post-departure epochs by basic probabilistic reasoning and discrete-time renewal theory. In steady-state, let ${\left\{P_j\right\}}_0^N$ and ${\left\{P_j^{-}\right\}}_0^N$ be the probabilities representing the number of jobs in the queue at random times and before arrival, respectively. Since the inter-arrival times use geometric distribution, the arrivals are independent of other events. Thus, it implies that $P_j=P_j^{-},\forall j=0,1,\dots ,N$; for details, see Boxma and Groenendijk [27]. If the server is idle, there are $<a$ jobs in the queue. Suppose $N_q$ is the number of tasks in the queue at some random time. At a random epoch, the steady-state probabilities are $P_{n,0}=P(N_q=n,\mathrm{server} \mathrm{idle}),0\le n\le a-1$, and $P_{n,1}=P(N_q=n,\mathrm{server} \mathrm{busy}),0\le n\le N$. Given this,

$$P_j=\left\{\begin{array}{cc}{P}_{j,0}+P_{j,1}\hfill & \mathrm{if} 0\le j\le a-1\hfill \\ P_{j,1}\hfill & \mathrm{if} a\le j\le N\hfill \end{array}\right.$$

Let the limiting pmf of the elapsed service time be ${\widehat{s}}_{\ell }$, which is determined by ${\widehat{s}}_{\ell }=\mu {\displaystyle \sum _{r=\ell +1}^{\infty }}{s}_r,\ell \ge 0$ (see, [5] (p. 20)), and ${\widehat{k}}_{\ell }$ be the probability that the number of arrivals during an elapsed service time is ℓ. This yields

$$\begin{array}{cc}\hfill {\widehat{k}}_{\ell }& ={\displaystyle \sum _{i=\ell }^{\infty }}\left(\genfrac{}{}{0pt}{}{i}{\ell }\right){\lambda }^{\ell }{(1-\lambda )}^{i-\ell }{\widehat{s}}_i,\ell =0,1,\dots \hfill \end{array}$$

If $E^{*}$ is the mean inter-departure time of processing batches, $1/E^{*}$ represents the departure rate. At the batch departure epoch, if the number of jobs in the queue is less than a, the subsequent batch departure occurs after an idle time and the service time processing time. Otherwise, the release of the next batch takes place following the processing time of the subsequent batch. This gives,

$$\begin{array}{cc}\hfill E^{*}& =E\left(S\right)\left(1-{\displaystyle \sum _{i=0}^{a-1}}{P}_i^{+}\right)+{\displaystyle \sum _{i=0}^{a-1}}{P}_i^{+}\left(\frac{a-i}{\lambda }+E\left(S\right)\right)=E\left(S\right)+{\displaystyle \sum _{i=0}^{a-1}}{P}_i^{+}\frac{(a-i)}{\lambda }.\hfill \end{array}$$

Remark 7.

It can also be put down as $E^{*}=E\left(S\right)+{\displaystyle \sum _{i=1}^a}{P}_{a-i}^{+}\left(\frac{i}{\lambda }\right)$. When $a=1$, it matches with Chaudhry and Chang [7].

Theorem 1.

The random- and post-departure-epoch probabilities $\{{\left\{P_{j,0}\right\}}_0^{a-1},{\left\{P_{j,1}\right\}}_0^N\}$, and ${\left\{P_j^{+}\right\}}_0^N$ are related by

$$\begin{array}{cc}\hfill P_{j,0}& =\frac{{\displaystyle \sum _{i=0}^j}{P}_i^{+}}{{\displaystyle \sum _{i=0}^{a-1}}(a-i)P_i^{+}+\lambda E\left(S\right)},0\le j\le a-1,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill P_{j,1}& =\left(1-{\displaystyle \sum _{i=0}^{a-1}}{P}_{i,0}\right)·[{\displaystyle \sum _{i=0}^{a-1}}{P}_i^{+}{\widehat{k}}_j+{\displaystyle \sum _{i=a}^{b-1}}{P}_i^{+}\left({\displaystyle \sum _{m=i}^b}{y}_m\right){\widehat{k}}_j+{\displaystyle \sum _{i=a+1}^{b-1}}{P}_i^{+}{\displaystyle \sum _{m=a}^{i-1}}{y}_m{\widehat{k}}_{j-i+m}\hfill \\ & +{\displaystyle \sum _{i=b}^N}{P}_i^{+}{\displaystyle \sum _{m=a}^b}{y}_m{\widehat{k}}_{j-i+m}],0\le j\le N-1.\hfill \end{array}$$

(12)

Finally, $P_{N,1}$ can be found from $P_{N,1}=1-{\displaystyle \sum _{j=0}^{a-1}}{P}_{j,0}-{\displaystyle \sum _{j=0}^{N-1}}{P}_{j,1}.$

Proof. 

The fraction of the time the batch server remains idle between two successive departure epochs is the probability of getting the server idle at a random epoch ($P_{idle}$). Let $E\left(I\right)$ be the mean idle period. Using the definition of $E\left(S\right)$ and $E\left(I\right)$, we have

$$\begin{array}{cc}\hfill P_{idle}& =\frac{E\left(I\right)}{E\left(I\right)+E\left(S\right)},\hfill \\ \hfill P_{j,0}& =P_{idle}\times P\left(\mathrm{fraction} \mathrm{of} \mathrm{idle} \mathrm{perio}\right)=\frac{E\left(I\right)}{E\left(I\right)+E\left(S\right)}\times \frac{\frac{1}{\lambda }{\displaystyle \sum _{i=0}^j}{P}_i^{+}}{E\left(I\right)},0\le j\le a-1,\hfill \end{array}$$

where $E\left(I\right)=\frac{1}{\lambda }{\displaystyle \sum _{i=0}^{a-1}}{P}_i^{+}(a-i)$. We employ system state conditioning and discreet renewal theory to find $P_{j,1}$. The processor is active with probability $(1-{\displaystyle \sum _{i=0}^{a-1}}{P}_{i,0})$; thus,

$$\begin{array}{cc}\hfill P_{j,1}& =P(N_q=j,\mathrm{processor} \mathrm{active})=(1-{\displaystyle \sum _{i=0}^{a-1}}{P}_{i,0})P(N_q=j|\mathrm{processor} \mathrm{active})\hfill \end{array}$$

(13)

Assuming that $k_j^e$ is the number of jobs that come following an embedded Markov point until the elapsed service time, we have

$$\begin{array}{cc}\hfill P(N_q=j|\mathrm{processor} \mathrm{activ})& ={\displaystyle \sum _{i=0}^{a-1}}{P}_i^{+}{\displaystyle \sum _{m=a}^b}{y}_m{\widehat{k}}_j+{\displaystyle \sum _{i=a}^{b-1}}{P}_i^{+}{\displaystyle \sum _{m=i}^b}{y}_m{\widehat{k}}_j\hfill \\ & +{\displaystyle \sum _{i=a+1}^{b-1}}{P}_i^{+}{\displaystyle \sum _{m=a}^{i-1}}{y}_m{\widehat{k}}_{j-i+m}+{\displaystyle \sum _{i=b}^N}{P}_i^{+}{\displaystyle \sum _{m=a}^b}{y}_m{\widehat{k}}_{j-i+m},\hfill \\ & 0\le j\le N-1.\hfill \end{array}$$

(14)

Putting together (13) and (14), we obtain (12). We can obtain $P_{N,1}$ using normalization condition. Thus, we obtain random epoch probabilities $\{{\left\{P_{j,0}\right\}}_0^{a-1},{\left\{P_{j,1}\right\}}_0^N\}$ in connection with post departure epoch probabilities ${\left\{P_j^{+}\right\}}_0^N$. □

Though the relations between random- and post-departure epoch probabilities are available in [13] using transition rates, here, we develop an alternative method to obtain the queue-length distributions at random epochs.

Because of BASTA (Bernoulli arrivals see time averages) property, see [5], the queue-length distribution exactly before arrival of job will be equal to that of $P_{j,0}$ and $P_{j,1}$. Further, since the outside observation epoch falls in an interval between arrival and departure epochs, the outside observer’s distribution is the same as the random epoch distribution.

4. Performance Measures

This section deals with several measures of performance. The average number of jobs in the queue ($L_q$) is given by

$$\begin{array}{cc}\hfill L_q& ={\displaystyle \sum _{j=0}^{a-1}}j{P}_{j,0}+{\displaystyle \sum _{j=0}^N}j{P}_{j,1}\hfill \end{array}$$

The probability of the processor being busy (PB) at some random moment is specified by $1-{\sum }_{j=0}^{a-1}{P}_{j,0}$. Due to the BASTA property, the loss or blocking probability ($PBL$) is given by $PBL=P_{N,1}$. Since the effective arrival rate ${\lambda }_e=\lambda (1-P_{N,1})$, we can obtain the average wait time in the queue ($W_q$) by employing Little’s law as $W_q=\frac{L_q}{{\lambda }_e}$. The reported result of $W_q$ in [13] is incorrect. They have applied the effective arrival rate as $\lambda$ instead of ${\lambda }_e$.

Remark 8.

If $\rho <1$ and $N\to \infty$, then ${\lambda }_e=\lambda$ and $L_q={\displaystyle \sum _{j=0}^{a-1}}j{P}_{j,0}+{\displaystyle \sum _{j=0}^{\infty }}j{P}_{j,1}$. Using Little’s law, the average waiting time in the queue ($W_q$) can be computed as $W_q=\frac{L_q}{\lambda }$.

5. Numerical Results

To exemplify the analytic results found in this article, we illustrate several numerical outcomes in tables and figures. We also give several performance measures, for instance, the average queue length ($L_q$), the average waiting time in the queue ($W_q$), and the probability of loss ($P_{loss}$). The computations were performed in double precision and executed in a 64-bit windows 10 professional OS possessing Intel(R) corei5-6200U processor $@2.30$ GHz and 8 GB DDR3 RAM utilizing MAPLE 18 software. The numerical results were exact, but we reported outcomes rounding to six decimal places. Because reported outcomes are rounded, the sum of probabilities may not add up to one in some cases.

Example 1.

$Geo/N{B}_2^{a,Y}/1/10$ queue. Consider the distribution of service time as being in two stage negative binomial distribution (NB) with PGF $S\left(z\right)={\left(\frac{\mu z}{1-\overline{\mu }z}\right)}^2$. Here, we consider the same parameters as in Table 1 of the paper [13] to compare the results. The parameters are $E\left(S\right)=5,$$y_2=0.2,y_3=0.1$, $y_4=0.7$, and $E\left(Y\right)=3.5$. The arrival rates are $0.14,0.35$, and $0.56$ with corresponding traffic intensities $\rho =0.2,0.5$, and $0.8$, respectively. To show the evaluation process, let us assume $\lambda =0.14$. The denominator of Equation (6) has six roots, two of which are outside, and four are in and on the unit circle. From Equation (8),

$$P_N^{+}\left(z\right)=T\left(0.068061+\frac{7.817631}{z-6.392624}-\frac{6.204397}{z-5.024572}\right),$$

where $T=6.250001$. Similarly, the denominator of $P_N^{+}\left(z\right)$ has two roots outside the unit circle when arrival rates are $0.35$ and $0.56$, and we have from Equation (8),

$$\begin{array}{c}\hfill P_N^{+}\left(z\right)=T\left(c_0+\frac{M_{b+1}}{z-{\gamma }_{b+1}}+\frac{M_{b+2}}{z-{\gamma }_{b+2}}\right),\end{array}$$

where ${\gamma }_{b+1}=2.106094$, ${\gamma }_{b+2}=3.454495$, $T=6.256038$, $c_0=0.053212,M_{b+1}=-0.619857$, $M_{b+2}=1.113392$ and ${\gamma }_{b+1}=1.304299$, ${\gamma }_{b+2}=2.688401$, $T=6.867372$, $c_0=0.025779$, $M_{b+1}=-0.081358$, $M_{b+2}=0.224798$, respectively.

Now, we can find post-departure epoch probabilities from Equation (9). The results are presented in

Table 1. Queue-length distributions at various epochs for the $Geo/N{B}_2^{2,Y}/1/10$ queue.

	${\mathit{P}}_{\mathit{j}}^{+}$				$\{{\mathit{P}}_{\mathit{j},0},{\mathit{P}}_{\mathit{j},1}\}$
	$\mathbf{\rho }=\mathbf{0}.\mathbf{2}$	$\mathbf{\rho }=\mathbf{0}.\mathbf{5}$	$\mathbf{\rho }=\mathbf{0}.\mathbf{8}$		$\mathbf{\rho }=\mathbf{0}.\mathbf{2}$	$\mathbf{\rho }=\mathbf{0}.\mathbf{5}$	$\mathbf{\rho }=\mathbf{0}.\mathbf{8}$
$P_0^{+}$	0.499734	0.157755	0.030728	$P_{0,0}$	0.244991	0.066960	0.010330
$P_1^{+}$	0.340335	0.290459	0.113199	$P_{1,0}$	0.411838	0.190246	0.048385
$P_2^{+}$	0.118658	0.246052	0.169892	$P_{0,1}$	0.239989	0.301616	0.170317
$P_3^{+}$	0.031582	0.148137	0.161190	$P_{1,1}$	0.077189	0.206250	0.172979
$P_4^{+}$	0.007532	0.079398	0.135131	$P_{2,1}$	0.019965	0.116185	0.148409
$P_5^{+}$	0.001694	0.040301	0.107666	$P_{3,1}$	0.004694	0.060454	0.119818
$P_6^{+}$	0.000368	0.019932	0.084302	$P_{4,1}$	0.001047	0.030223	0.093880
$P_7^{+}$	0.000078	0.009692	0.065970	$P_{5,1}$	0.000226	0.014762	0.071726
$P_8^{+}$	0.000016	0.004598	0.049310	$P_{6,1}$	0.000048	0.007247	0.059329
$P_9^{+}$	0.000003	0.002109	0.033773	$P_{7,1}$	0.000010	0.003379	0.041064
$P_{10}^{+}$	0.000001	0.001567	0.048838	$P_{8,1}$	0.000002	0.001560	0.027772
				$P_{9,1}$	0.000000	0.000665	0.016418
				$P_{10,1}$	0.000000	0.000452	0.019571
Sum	1.000000	1.000000	1.000000	Sum	1.000000	1.000000	1.000000
				$L_q$	0.548733	1.095056	2.820867
				$W_q$	3.919521	3.130145	5.137817
				$PBL$	0.000000	0.000452	0.019571

. We note that the results of queue-length distribution at post-departure epoch match the results given by Yi et al. [13], but the random epoch does not. We have also computed queue-length distributions at random epochs using their method for checking purposes. They match perfectly with our results. However, the results presented in the paper by [13] are different. Thus, various performance measures are also not the same.

Remark 9.

It may be noted that $P_N^{+}\left(z\right)$ is a polynomial in both cases, as can be seen in Table 1. The same applies to other cases as well.

Example 2.

$Geo/DP{H}^{a,Y}/1/N$ queue. The service-time distribution is assumed to be discrete phase-type (DPH) having $s_k=\mathit{\alpha }{\mathit{T}}^{\mathit{k}-\mathbf{1}}{\mathit{T}}^{\mathbf{0}}$, $k=1,2,\dots$, ${\mathit{T}}^{\mathbf{0}}=\mathit{e}-T\mathit{e}$, where $\mathit{e}$ is the appropriate column vector with all elements equal in size. This gives the PGF of service-time distribution as $S\left(z\right)=\mathit{z}\mathit{\alpha }{(\mathit{I}-\mathit{zT})}^{-\mathbf{1}}{\mathit{T}}^{\mathbf{0}}$, $\left|z\right|\le 1$.

Table 2. Queue-length distributions at various epochs for the $Geo/DP{H}^{3,Y}/1/15$ and $Geo/DP{H}^{3,Y}/1/50$ queues.

	$Geo/DP{H}^{3,Y}/1/15$				$Geo/DP{H}^{3,Y}/1/50$
	${\mathit{y}}_{\mathbf{3}}=\mathbf{0}.\mathbf{5},{\mathit{y}}_{\mathbf{5}}=\mathbf{0}.\mathbf{3},{\mathit{y}}_{\mathbf{8}}=\mathbf{0}.\mathbf{2}$				${\mathit{y}}_{\mathbf{3}}=\mathbf{0}.\mathbf{7},{\mathit{y}}_{\mathbf{4}}=\mathbf{0}.\mathbf{2},{\mathit{y}}_{\mathbf{5}}=\mathbf{0}.\mathbf{1}$
	$\mathit{\lambda }=\mathbf{0}.\mathbf{7},\mathit{\rho }=\mathbf{0}.\mathbf{30515}$				$\mathit{\lambda }=\mathbf{0}.\mathbf{84},\mathit{\rho }=\mathbf{1}.\mathbf{0376}$
$\mathit{j}$	${\mathit{P}}_{\mathit{j}}^{+}$	${\mathit{P}}_{\mathit{j},\mathbf{0}}$	${\mathit{P}}_{\mathit{j},\mathbf{1}}$	$\mathit{j}$	${\mathit{P}}_{\mathit{j}}^{+}$	${\mathit{P}}_{\mathit{j},\mathbf{0}}$	${\mathit{P}}_{\mathit{j},\mathbf{1}}$
0	0.170394	0.055766	0.258523	0	0.001537	0.000432	0.009680
1	0.469965	0.209576	0.111943	1	0.009682	0.00315	0.009988
2	0.200698	0.27526	0.049464	2	0.010114	0.005989	0.010254
3	0.088375		0.021963	3	0.010392		0.010524
4	0.039275		0.00974	4	0.010666		0.010801
5	0.017414		0.004321	5	0.010946		0.011085
10	0.000299		0.000074	10	0.012463		0.012621
11	0.000133		0.000032	20	0.016158		0.016363
12	0.000059		0.000014	30	0.020947		0.021213
13	0.000026		0.000005	40	0.027156		0.027501
14	0.000012		0.000001	49	0.034304		0.009883
15	0.000005		0.000004	50	0.035206		0.049366
Sum	1.000000	0.540603	0.459397	Sum	1.000000	0.009570	0.990430
	$L_q$ = 1.1208, $W_q$ = 1.6012				$L_q$ = 29.6650, $W_q$ = 7.1493
	$PBL$ = 0.000004				$PBL$ = 0.049366

shows the queue length distributions at different times employing the DPH service time distribution. For the first example of Table 2, we suppose

$$\mathit{\alpha }=\left[\begin{array}{ccc}0.40& 0.50& 0.10\end{array}\right], \mathbf{T}=\left[\begin{array}{ccc}0.10& 0.20& 0.05\\ 0.30& 0.15& 0.10\\ 0.20& 0.50& 0.10\end{array}\right],$$

with $E\left(S\right)=2.005267$ and the other parameters are $\lambda =0.7,N=15,y_3=0.5,y_5=0.3,$ $y_8=0.2$ with $E\left(Y\right)=4.6$. Here the denominator of $P_N^{+}\left(z\right)$ has eleven distinct roots, out of which three roots are outside the unit circle, and they are ${\gamma }_{b+1}=2.254374$, ${\gamma }_{b+2}=-9.788249$, and ${\gamma }_{b+3}=-222.447161$. From Equation (8),

$$\begin{array}{cc}\hfill P_N^{+}\left(z\right)& =T\left(-3.63007+\frac{1.301912}{z-2.254374}-\frac{0.702111}{z+9.788249}+\frac{930.333318}{z+222.447161}\right),\hfill \end{array}$$

(15)

where $T=-1.755933$.

In the second example of Table 2, we assume $\mathit{\alpha }=\left[\begin{array}{cc}0.60& 0.40\end{array}\right]$, and $\mathit{T}=\left[\begin{array}{cc}0.5& 1/3\\ 1/3& 1/3\end{array}\right]$ with $E\left(S\right)=4.2$ and the other parameters are $\lambda =0.84,N=50,y_3=0.7,y_4=0.2,y_5=0.1$ with $E\left(Y\right)=3.4$. As $\rho =1.0376$, we have one root in the range $(0,1)$ from the remaining four roots, and the other three are outside the circle of unity. From Equation (8),

$$\begin{array}{cc}\hfill P_N^{+}\left(z\right)& =T\left(0.000674+\frac{5.2\times {10}^{-21}}{z-1.375407}+\frac{0.001238}{z-0.974373}-\frac{0.006324}{z-16.101113}-\frac{4.900051\times {10}^{-13}}{z-16.100784}\right),\hfill \end{array}$$

(16)

where $T=-7.569441$.

Remark 10.

It may be noted that Equations (15) and (16) are polynomials since the coefficients of $z^b$ are zero. This may be seen in Table 1. This also applies to all the examples that follow.

Example 3.

As before, in this example, we consider two cases: (i) $Geo/MGe{o}_2^{2,Y}/1/10$ queue. Here, the service time is a mixture of two geometric distributions with PGF $S\left(z\right)={\varsigma }_1\frac{{\mu }_{1}z}{1-(1-{\mu }_1)z}+{\varsigma }_2\frac{{\mu }_{2}z}{1-(1-{\mu }_2)z}$, where ${\varsigma }_1+{\varsigma }_2=1$. The parameters are taken as $E\left(S\right)=2.83333,N=10,{\varsigma }_1=0.6,{\mu }_1=0.4,{\mu }_2=0.3,y_2=0.2,y_3=0.1,y_4=0.7,E\left(Y\right)=3.5$, $\lambda =0.14$. The denominator of $P_N^{+}\left(z\right)$ has six distinct roots, out of which two roots are outside the unit circle, and they are ${\gamma }_{b+1}=4.031316$ and ${\gamma }_{b+2}=5.727647$. From Equation (8),

$$\begin{array}{cc}\hfill P_N^{+}\left(z\right)& =T\left(-0.561775-\frac{1.819848}{z-4.031316}-\frac{4.545280}{z-5.727647}\right),\hfill \end{array}$$

where $T=1.0000001$.

(ii) $Geo/D^{2,Y}/1/50$ queue. Here, we consider service-time distribution as deterministic with PGF $S\left(z\right)=z^k$ for some constant $k=4,s_k=1$. The parameters are taken as $E\left(S\right)=4,N=50,y_2=0.5,y_3=0.2,y_4=0.1,y_5=0.1,y_6=0.1,E\left(Y\right)=3.1$, and $\lambda =0.4$. In this case, the denominator of $P_N^{+}\left(z\right)$ has eight distinct roots, out of which two roots are outside the unit circle, and they are ${\gamma }_{b+1}=5.019701$ and ${\gamma }_{b+2}=-11.795711$. From Equation (8),

$$\begin{array}{cc}\hfill P_N^{+}\left(z\right)& =T\left(-1.316734z^2+1.021799z-102.664691+\frac{879.864750}{z+11.795711}-\frac{141.482007}{z-5.019701}\right),\hfill \end{array}$$

where $T=1$.

Table 3. Queue-length distributions at various epochs for the $Geo/MGe{o}_2^{2,Y}/1/10$ and $Geo/D^{2,Y}/1/50$ queues.

	$\mathit{Geo}/{\mathit{MGeo}}_2^{2,\mathit{Y}}/1/10$				$\mathit{Geo}/{\mathit{D}}^{2,\mathit{Y}}/1/50$
	${\mathit{y}}_{\mathbf{2}}=\mathbf{0}.\mathbf{2},{\mathit{y}}_{\mathbf{3}}=\mathbf{0}.\mathbf{1},{\mathit{y}}_{\mathbf{4}}=\mathbf{0}.\mathbf{7}$				${\mathit{y}}_{\mathbf{2}}=\mathbf{0}.\mathbf{5},{\mathit{y}}_{\mathbf{3}}=\mathbf{0}.\mathbf{2},{\mathit{y}}_{\mathbf{4}}=\mathbf{0}.\mathbf{1}$
	${\mathit{\mu }}_{\mathbf{1}}=\mathbf{0}.\mathbf{4},{\mathit{\mu }}_{\mathbf{2}}=\mathbf{0}.\mathbf{3}$,${\mathit{\varsigma }}_{\mathbf{1}}=\mathbf{0}.\mathbf{6}$				${\mathit{y}}_{\mathbf{5}}=\mathbf{0}.\mathbf{1},{\mathit{y}}_{\mathbf{6}}=\mathbf{0}.\mathbf{1}$
	$\mathit{\lambda }=\mathbf{0}.\mathbf{14},\mathit{\rho }=\mathbf{0}.\mathbf{11333}$				$\mathit{\lambda }=\mathbf{0}.\mathbf{4},\mathit{\rho }=\mathbf{0}.\mathbf{51613}$
$\mathit{j}$	${\mathit{P}}_{\mathit{j}}^{+}$	${\mathit{P}}_{\mathit{j},\mathbf{0}}$	${\mathit{P}}_{\mathit{j},\mathbf{1}}$	$\mathit{j}$	${\mathit{P}}_{\mathit{j}}^{+}$	${\mathit{P}}_{\mathit{j},\mathbf{0}}$	${\mathit{P}}_{\mathit{j},\mathbf{1}}$
0	0.683222	0.339297	0.15549	0	0.112580	0.052650	0.353602
1	0.250531	0.463714	0.032491	1	0.313097	0.199077	0.253667
2	0.051967		0.006998	2	0.337945		0.107802
3	0.011114		0.001551	3	0.177390		0.027021
4	0.002447		0.000352	4	0.048246		0.004913
5	0.000553		0.000082	5	0.008517		0.001019
6	0.000128		0.000019	10	0.000003		0.000003
7	0.000030		0.000005	20	0.000000		0.000000
8	0.000007		0.000001	30	0.000000		0.000000
9	0.000002		0.000000	40	0.000000		0.000000
10	0.000001		0.000000	50	0.000000		0.000000
Sum	1.000000	0.803010	0.196990	Sum	1.000000	0.251727	0.748273
	$L_q$ = 0.516831, $W_q$ = 3.691651				$L_q$ = 0.775713, $W_q$ = 1.939284
	$PBL$ = 0.000000				$PBL$ = 0.000000

presents queue-length distributions at various epochs when the service-time distributions are a mixture of geometric and deterministic.

Example 4.

Here, we consider two cases.

Table 4. Queue-length distributions at various epochs when $N\to \infty$.

	$\mathit{Geo}/{\mathit{NB}}^{2,\mathit{Y}}/1/\mathit{\infty }$			$\mathit{Geo}/{\mathit{DPH}}^{2,\mathit{Y}}/1/\mathit{\infty }$
	${\mathit{y}}_{\mathbf{2}}=\mathbf{0}.\mathbf{5},{\mathit{y}}_{\mathbf{4}}=\mathbf{0}.\mathbf{2},{\mathit{y}}_{\mathbf{5}}=\mathbf{0}.\mathbf{1},{\mathit{y}}_{\mathbf{6}}=\mathbf{0}.\mathbf{1}$			${\mathit{y}}_{\mathbf{2}}=\mathbf{0}.\mathbf{5},{\mathit{y}}_{\mathbf{3}}=\mathbf{0}.\mathbf{3},{\mathit{y}}_{\mathbf{4}}=\mathbf{0}.\mathbf{2},$
	${\mathit{y}}_{\mathbf{8}}=\mathbf{0}.\mathbf{1},\mathit{\lambda }=\mathbf{0}.\mathbf{703},\mathit{\rho }=\mathbf{0}.\mathbf{95}$			$\mathit{\lambda }=\mathbf{0}.\mathbf{58},\mathit{\rho }=\mathbf{0}.\mathbf{902}$
$\mathit{j}$	${\mathit{P}}_{\mathit{j}}^{+}$	${\mathit{P}}_{\mathit{j},\mathbf{0}}$	${\mathit{P}}_{\mathit{j},\mathbf{1}}$	${\mathit{P}}_{\mathit{j}}^{+}$	${\mathit{P}}_{\mathit{j},\mathbf{0}}$	${\mathit{P}}_{\mathit{j},\mathbf{1}}$
0	0.002177	0.000616	0.028874	0.030809	0.011992	0.070392
1	0.013258	0.004369	0.034832	0.071536	0.039836	0.065422
2	0.029345		0.036239	0.066885		0.060557
3	0.034741		0.035880	0.061929		0.056044
4	0.035937		0.034858	0.057314		0.051866
5	0.035505		0.033609	0.053041		0.047999
10	0.029400		0.027374	0.036007		0.032584
20	0.019366		0.018028	0.016594		0.015016
30	0.012753		0.011872	0.007647		0.006920
50	0.005531		0.005149	0.001624		0.001470
100	0.000685		0.000638	0.000034		0.000031
150	0.000085		0.000079	0.000001		0.000001
200	0.000011		0.000010	0.000000		0.000000
250	0.000001		0.000001	0.000000		0.000000
300	0.000000		0.000000	0.000000		0.000000
500	0.000000		0.000000	0.000000		0.000000
Sum	1.000000	0.004986	0.995014	1.000000	0.051828	0.948172
	$L_q$ =23.80989, $W_q$ =33.86898			$L_q$ =11.815062, $W_q$ =20.370796

presents the results of $Geo/G^{a,Y}/1/\infty$ which can be found from the $Geo/G^{a,Y}/1/N$ system by assuming $\rho <1$ and buffer capacity N appropriately large. We can easily compute the queue-length distributions of infinite queue capacity from finite queue capacity by assuming $\rho <1$ and N sufficiently large.

(i) $Geo/N{B}^{a,Y}/1/\infty$ queue. Here, we assume negative binomial (NB) service time distribution and the parameters are taken as $\lambda =0.703,E\left(S\right)=5,y_2=0.5,y_4=0.2,y_5=0.1,y_6=0.1,y_8=0.1$ with $E\left(Y\right)=3.7$. From Equation (8), we have

$$\begin{array}{cc}\hfill P_N^{+}\left(z\right)& =T\left(0.005129+\frac{0.03006}{z-2.520637}-\frac{0.00745}{z-1.042657}\right),\hfill \end{array}$$

where $T=6.250024$.

(ii) $Geo/DP{H}^{a,Y}/1/\infty$ queue. For a DPH service time distribution, the settings are chosenas

$$\mathit{\alpha }=\left[\begin{array}{cc}0.60& 0.40\end{array}\right] and \mathbf{T}=\left[\begin{array}{cc}0.5& 1/3\\ 1/3& 1/3\end{array}\right].$$

with $E\left(S\right)=4.2$ and the other parameters are $\lambda =0.58,y_2=0.5,y_3=0.3,y_4=0.2$, and $E\left(Y\right)=2.7$. From Equation (8), we have

$$\begin{array}{cc}\hfill P_N^{+}\left(z\right)& =T\left(-0.001463+\frac{0.019982}{z-22.872839}-\frac{3.590919\times {10}^{-14}}{z-22.870101}-\frac{0.004169}{z-1.080550}+\frac{1.09\times {10}^{-20}}{z-1.543692}\right),\hfill \end{array}$$

where $T=20.250$.

In Figure 2, we compare the processing times to calculate probabilities at post-departure using the proposed technique and the method used by [13] (solving a linear system of equations (SLSE)) against finite buffer capacity. We take the NB service-time distribution in two stages with the input parameters in the same way as in Table 1 for $\rho =0.2$. We notice that, with the increase of N, the time needed by the roots method remains almost static, whereas the method used by [13] takes more time for larger buffer size N and increases remarkably as N increases. In the roots method, the variation of the time is minimal. This is because of the initial estimation of the Newton iterative method for the calculation of polynomial roots. The application of the SLSE process is low in reliability and time-consuming. The root method provides a faster solution and superior performance in solving a linear system of method equations both in speed and reliability.

Figure 3 shows the roots of the characteristic equation for the number in the queue with NB service time distribution having 4 successes is the convolution of 4 geometric distributions. Here we consider the parameters as $\lambda =0.81$, $y_4=0.4,y_6=0.1,y_{12}=0.1,y_{25}=0.2,y_{30}=0.1,y_{36}=0.1$, and $\rho =0.54$. There are 40 of roots inside, on, and out of the unit circle for the assumed parameters. Here the characteristic equation is

$$\begin{array}{cc}& z^{36}{\left(-0.886+0.486z\right)}^4-0.0256{\left(0.19+0.81z\right)}^4\left(0.1+0.4z^{32}+0.1z^{30}\right.\\ & \left.+0.1z^{24}+0.2z^{11}+0.1z^6\right).\end{array}$$

(17)

6. Conclusions

This article focuses on the $Geo/G^{a,Y}/1/N$ queue length distributions at various points in time. We use the roots of the associated characteristic equation to determine a unified way to compute performance measures for both infinite- and finite-buffer systems. Queue length distributions at a post-departure time are computed using an embedded Markov chain method. We obtain associations between queue length distributions at several time points by applying system state conditioning and discrete renewal theory. Several performance indices have been carried out, such as the blocking probability, the average wait time in the queue, and the average number of jobs in the queue. We illustrate them by using different numerical outcomes. The approach discussed in this paper can be used to cover cases when customers arrive in groups or even when arrivals are correlated (D-MAP-discrete Markovian arrival process).