**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–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–13]. In Gupta and Goswami [10], they discuss a discrete-time finitebuffer 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/*Ga*,*<sup>Y</sup>*/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/*Ga*,*<sup>b</sup>*/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 matrixgeometric method. In this connection, see the remark below. Following Chaudhry et al. [16],

**Citation:** Chaudhry, M.; Goswami, V. The *Geo*/*Ga*,*<sup>Y</sup>*/1/*N* Queue Revisited. *Mathematics* **2022**, *10*, 3142. https:// doi.org/10.3390/math10173142

Academic Editors: Gurami Tsitsiashvili and Alexander Bochkov

Received: 25 July 2022 Accepted: 23 August 2022 Published: 1 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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*/*Dc*/1, states that when *c* ≥ 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*/*Ga*,*<sup>Y</sup>*/1/*<sup>N</sup>* 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:


The remaining paper is structured as follows. Section 2 specifies the model. Section 3 analyzes the *Geo*/*Ga*,*<sup>Y</sup>*/1/*<sup>N</sup>* 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*/*Ga*,*<sup>Y</sup>*/1/*<sup>N</sup>* queue where jobs arrive following a Bernoulli process with parameter *λ*. 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*) = *yi*, *i* = *a*, *a* + 1, ... , *b* with ∑*bi*=*<sup>a</sup> yi* = 1, the probability generating function (PGF) *<sup>Y</sup>*(*z*) = ∑*bi*=*<sup>a</sup> yizi*, mean *E*(*Y*) = ∑*bi*=*<sup>a</sup> iyi* and *Y*-(*z*) = *<sup>z</sup>bY*(*z*<sup>−</sup><sup>1</sup>). 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 *yi* (when there are *Y* ≤ *b*, it takes all of them). When the number of jobs comes down below a threshold value *a*(≥ 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 {*Sn*, *n* ≥ 1} are independently and identically distributed (iid) with arbitrary pmf *sk* = *P*(*Sn* = *k*), *k* = 1, 2, ... and *s*0 = 0, the PGF *<sup>S</sup>*(*z*) = ∑∞*<sup>i</sup>*=<sup>1</sup> *siz<sup>i</sup>* and mean service time *E*(*S*) = *s* = *ddz<sup>S</sup>*(*z*)|*<sup>z</sup>*=1<sup>=</sup> 1/*μ*.

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* ≤ *N*. Thus, in the system, no more than *N* + *b* jobs can be available anytime. We presume offered load of the system as *ρ* = *<sup>λ</sup>E*(*S*)/*E*(*Y*). In LAS-DA policy, arrivals occur in (*<sup>u</sup>*<sup>−</sup>, *<sup>u</sup>*), and departures take place in (*<sup>u</sup>*, *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].

(u+,(<sup>u</sup> +

:Potential arrival epoch :Potential departure epoch

1)−):Outside observer's interval u<sup>+</sup> :Epoch after a potential departure

∗:Outside observer's epoch u− :Epochprior toapotentialarrival

**Figure 1.** Various time epochs in LAS-DA.

## **3. Queue-Length Distributions**

Here, we find steady-state queue-length distributions at various epochs of *Geo*/*Ga*,*<sup>Y</sup>*/1/*<sup>N</sup>* queue with the LAS-DA policy.

#### *3.1. Post-Departure Epoch Probabilities*

Let *Q*<sup>+</sup>- be the number of jobs in the queue after completing the th service. Suppose *A*-+1 and *Y*-+1 represent the number of arrivals throughout the processing time on the (- + 1)th job and the processing capacity of (- + 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)\_{\prime},$$

where *<sup>x</sup>*+=max(*<sup>x</sup>*, <sup>0</sup>). The probability distribution of *A*-+1 is

$$\begin{aligned} k\_n &= P(A\_{\ell+1} = n) = \sum\_{j=1}^{\infty} P(A\_{\ell+1} = n | S\_{\ell+1} = j) P(S\_{\ell+1} = j), \\ &= \sum\_{j=n}^{\infty} s\_j \binom{j}{n} \lambda^n (1 - \lambda)^{j-n}, \ n \ge 0. \end{aligned}$$

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 (- + 1)th job is *j* slots, then there will be *j* time slots where *n* arrivals may occur. One may note that (*j n*)*λ<sup>n</sup>*(<sup>1</sup> − *<sup>λ</sup>*)*j*−*<sup>n</sup>* is the probability of *n* arrivals in *j* slots. Let *<sup>K</sup>*(*z*) = ∑ ∞ *<sup>n</sup>*=0 *knzn* be the probability generating function of the sequence {*kn*, *n* = 0, 1, . . .}. Thus

$$K(z) = \sum\_{n=0}^{\infty} \sum\_{j=n}^{\infty} s\_j \binom{j}{n} \lambda^n (1-\lambda)^{j-n} z^n = \sum\_{j=0}^{\infty} s\_j (1-\lambda+\lambda z)^j = S(1-\lambda+\lambda z).$$

Transition probabilities in one step of underlying Markov chain for *pij* = lim -→∞*Pr*{*Q*<sup>+</sup> -+1 = *j*|*Q*<sup>+</sup> - = *i*} are given as

$$p\_{ij} = \begin{cases} k\_{jr} & 0 \le i \le a, \ 0 \le j \le N - 1, \\ k\_j \sum\_{r=i}^b y\_r + \sum\_{r=a}^{i-1} y\_r k\_{j-i+r}, & a+1 \le i \le b-1, \ i-a-1 \le j \le N-1, \\ \sum\_{r=a}^b y\_r k\_{j-i+r}, & b \le i \le N, \ i-b \le j \le N-1, \\ \ell\_{N\_r} & 0 \le i \le a, \ j = N\_r \\ \ell\_N \sum\_{r=i}^b y\_r + \sum\_{r=a}^{i-1} y\_r \ell\_{N-i+r}, & a+1 \le i \le b-1, \ j = N\_r \\ \sum\_{r=a}^b y\_r \ell\_{N-i+r} & b \le i \le N, \ j = N, \end{cases} \tag{1}$$

where *j* = ∞ ∑ *<sup>r</sup>*=*j kr* and *kj* with *j* < 0 defined to be zero, and which leads to the transition probability matrix *P* = (*pij*) as

*P* = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ 0 1 ... *j* ... *N* − 1 *N* 0 *k*0 *k*1 ... *kj* ... *kN*−<sup>1</sup> -*N* 1 *k*0 *k*1 ... *kj* ... *kN*−<sup>1</sup> -*N* . . . . . . . . . . . . . . . . . . . . . . . . *a k*0 *k*1 . . . *kj* ... *kN*−<sup>1</sup> -*N a* + 1 *k*0 *b* ∑ *<sup>r</sup>*=*a*+1 *yr k*1 *b* ∑ *<sup>r</sup>*=*a*+1 *yr* . . . *kj b* ∑ *<sup>r</sup>*=*a*+1 *yr* ... *kN*−<sup>1</sup> *b* ∑ *<sup>r</sup>*=*a*+1 *yr* -*N b* ∑ *<sup>r</sup>*=*a*+1 *yr* <sup>+</sup>*k*0*ya* <sup>+</sup>*kj*−<sup>1</sup>*ya* <sup>+</sup>*kN*−<sup>2</sup>*ya* +-*<sup>N</sup>*−<sup>1</sup>*ya* . . . . . . . . . . . . . . . . . . . . . . . . *b k*0*yb b* ∑ *<sup>r</sup>*=*b*−1 *yrk*1−*b*+*<sup>r</sup>* . . . *b* ∑ *r*=*a yrkj*−*b*+*<sup>r</sup>* ... *b* ∑ *r*=*a yrkN*−1−*b*+*<sup>r</sup> b* ∑ *r*=*a yrkN*−*b*+*<sup>r</sup> b* + 1 0 *k*0*yb* . . . *b* ∑ *r*=*a yrkj*−*b*+*r*−<sup>1</sup> ... *b* ∑ *r*=*a yrkN*−*b*+*r*−<sup>1</sup> *b* ∑ *r*=*a yrkN*−*b*+*r*−<sup>1</sup> . . . . . . . . . . . . . . . . . . . . . . . . *N* 0 00 *b* ∑ *r*=*a yrkj*−*N*+*<sup>r</sup>* . . . *b* ∑ *r*=*a yrkr*−<sup>1</sup> *b* ∑ *r*=*a yrr* ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠

.

The (- + 1)th batch service starts when there are '*a*' jobs in the queue, and the state changeover occurs from 0 ≤ *i* ≤ *a* − 1 to 0 ≤ *j* ≤ *N* − 1. If the state changeover is from state *i* ≥ *a* to state 0 ≤ *j* ≤ *N* − 1, then there is a busy period between the leaving epoch of the th batch and the commencement of processing (- + <sup>1</sup>)st batch. Let the steady-state probability *p*+ = {*P*+0 , *P*+1 , ... , *<sup>P</sup>*+*n* , ... , *P*+*N* } represent - jobs at departure epochs. Then, in steady-state, *p*+ = *p*<sup>+</sup>*<sup>P</sup>* can be expressed as follows:

$$\begin{aligned} P\_j^+ &= \sum\_{i=0}^a P\_i^+ k\_j + \sum\_{i=a+1}^{b-1} P\_i^+ \sum\_{r=i}^b y\_r k\_j + \sum\_{i=a+1}^{b-1} P\_i^+ \sum\_{r=a}^{i-1} y\_r k\_{j-i+r} + \sum\_{i=b}^N P\_i^+ \sum\_{r=a}^b y\_r k\_{j-i+r} \\ &0 \le j \le N - 1, \end{aligned} \tag{2}$$

$$P\_N^+ = \sum\_{i=0}^a P\_i^+ \ell\_N + \sum\_{i=a+1}^{b-1} P\_i^+ \sum\_{r=i}^b y\_r \ell\_N + \sum\_{i=a+1}^{b-1} P\_i^+ \sum\_{r=a}^{i-1} y\_r \ell\_{N-i+r} + \sum\_{i=b}^N P\_i^+ \sum\_{r=a}^b y\_r \ell\_{N-i+r}, \tag{3}$$

where the normalization condition is *N* ∑ *j*=0 *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 *<sup>P</sup>*+(*z*) = *N* ∑ *j*=0 *P*+*j <sup>z</sup>j*. Multiplying Equation (2) by *zj* and then adding overall *j*, we obtain

$$\begin{split} P^{+}(z) &= \sum\_{j=0}^{N-1} k\_{j} z^{j} \sum\_{i=0}^{a} P^{+}\_{i} + \sum\_{j=0}^{N-1} z^{j} \sum\_{i=a+1}^{b-1} P^{+}\_{i} \left( \sum\_{r=i}^{b} y\_{r} \right) k\_{j} \\ &+ \sum\_{j=0}^{N-1} z^{j} \sum\_{i=a+1}^{b-1} P^{+}\_{i} \sum\_{r=a}^{i-1} y\_{r} k\_{j-i+r} + \sum\_{j=0}^{N-1} z^{j} \sum\_{i=b}^{N} P^{+}\_{i} \sum\_{r=a}^{b} y\_{r} k\_{j-i+r} + P^{+}\_{N} z^{N} \\ &\quad P^{+}(z) \left[ 1 - K(z) Y \left( \frac{1}{z} \right) \right] = K(z) \sum\_{i=0}^{a-1} P^{+}\_{i} \left( 1 - z^{i} Y \left( \frac{1}{z} \right) \right) + P^{+}\_{N} z^{N} \\ &\quad + K(z) \sum\_{i=a}^{b} P^{+}\_{i} \left( \sum\_{r=i}^{b} y\_{r} - z^{i} \sum\_{r=i}^{b} z^{r} \right) - \sum\_{i=a}^{b} y\_{i} \sum\_{j=i+1}^{N} P^{+}\_{j} + \sum\_{r=N-j+a}^{\infty} k\_{r} z^{r+j-i} \\ &\quad - \left( \sum\_{i=0}^{a-1} p\_{i}^{+} + \sum\_{i=a}^{b} P^{+}\_{i} \sum\_{r=i}^{b} y\_{r} \right) \sum\_{j=N}^{\infty} k\_{j} z^{j} .\end{split}$$

Simplifying the above equation, we ge<sup>t</sup>

$$\begin{split} P^{+}(z) &= \frac{K(z)\left[\sum\_{i=0}^{a-1} P\_{i}^{+}\left(z^{b} - z^{i}Y'(z)\right) + \sum\_{i=a}^{b-1} P\_{i}^{+}\left(\sum\_{r=i}^{b} y\_{7}z^{b} - z^{i}\sum\_{r=i}^{b} y\_{7}z^{b-r}\right)\right]}{z^{b} - K(z)Y'(z)} \\ &+ \frac{z^{N+b}\left(P\_{N}^{+} - \sum\_{i=a}^{b} y\_{i}\sum\_{j=i+1}^{N} P\_{j}^{+}\sum\_{r=N-j+a}^{\infty} k\_{7}z^{r+j-i-N}\right)}{z^{b} - K(z)Y'(z)} \\ &- \frac{z^{N+b}\left(\sum\_{i=0}^{a-1} P\_{i}^{+} + \sum\_{i=a}^{b} P\_{i}^{+}\sum\_{r=i}^{b} y\_{r}\right)\sum\_{j=N}^{\infty} k\_{j}z^{j-N}}{z^{b} - K(z)Y'(z)}. \end{split} \tag{4}$$

Only the first expression on the right side of the Equation (4) will add to the coefficient of *zj*, *j* = 0, 1, ... *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 *zj* for *j* ≤ *N* on both sides in Equation (4) to find *P*+*j* for *j* = 0, 1, . . . , *N*. Let

$$P\_N^+(z) = \frac{K(z)\left[\sum\_{i=0}^{a-1} P\_i^+\left(z^b - z^i Y'(z)\right) + \sum\_{i=a}^{b-1} P\_i^+\left(\sum\_{r=i}^b y\_r z^b - z^i \sum\_{r=i}^b y\_r z^{b-r}\right)\right]}{z^b - K(z)Y'(z)},\tag{5}$$

which is equivalent to the PGF of an infinite buffer case. The function *P*+*N* (*z*) is fully determined once *P*+*i* , *i* = 0, 1, ... , *b* − 1 are known. One may observe that when *ρ* < 1, Equation (5) denotes the PGF of discrete-time infinite buffer *Geo*/*Ga*,*<sup>Y</sup>*/1 queue. We can calculate the probabilities for *ρ* ≥ 1 in the case of a finite buffer *Geo*/*Ga*,*<sup>Y</sup>*/1/*<sup>N</sup>* queue.

**Remark 2.** *Using a* = 1*, the model is reduced to Geo*/*G<sup>Y</sup>*/1/*N queue.*

**Remark 3.** *Taking y*1 = 1*, yi* = 0, ∀ 2 ≤ *i* ≤ *b and a* = 1*, the model becomes Geo*/*G*/1/*N queue and Equation* (5) *establishes PGF as P*+*N* (*z*) = *P*+0 *<sup>K</sup>*(*z*)(<sup>1</sup>−*<sup>z</sup>*) *<sup>K</sup>*(*z*)−*<sup>z</sup> , where <sup>K</sup>*(*z*) = *S*(1 − *λ* + *<sup>λ</sup>z*)*, which corresponds to the results of [23].*

**Remark 4.** *Taking yb* = 1 *and yi* = 0, ∀ *i* = *b, the model is reduced to Geo*/*G*(*<sup>a</sup>*,*b*)/1/*<sup>N</sup> and Equation* (5) *establishes the PGF as P*+*N* (*z*) = *<sup>K</sup>*(*z*) *b*−1 ∑ *i*=0 *P*+*i* (*zb*−*z<sup>i</sup>*) *<sup>z</sup>b*−*<sup>K</sup>*(*z*) *.*

Intending to establish a unified method to solve the queueing system *Geo*/*Ga*,*<sup>Y</sup>*/1/*N*, we obtain {*P*+*n* }*N*0 from *<sup>P</sup>*+(*z*) by using the roots of characteristic equation and partialfraction 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 *<sup>K</sup>*(*z*) = *S*(1 − *λ* + *<sup>λ</sup>z*) as a rational function in *z*, specified by

$$K(z) = S(1 - \lambda + \lambda z) = \frac{f(z)}{g(z)},$$

where *f*(*z*) and *g*(*z*) 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),

$$P\_N^+(z) = \frac{f(z)\left[\sum\_{i=0}^{a-1} P\_i^+\left(z^b - z^i Y'(z)\right) + \sum\_{i=a}^{b-1} P\_i^+\left(\sum\_{r=i}^b y\_r z^b - z^i \sum\_{r=i}^b y\_r z^{b-r}\right)\right]}{z^b g(z) - f(z)Y'(z)}.\tag{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 | *z* |= 1, say, *γ*1 = 1, *γi* (*i* = 2, 3, . . . , *b* + *s*).

**Remark 5.** *If the denominator zbg*(*z*) − *f*(*z*)*Y*-(*z*) *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*

$$
\mu(y) = (y-2)(y-5)^2(y-7)^3(y-11).
$$

*restart* : *Digits* := 10 : *with*(*RootFinding*) :

*m* := (*y* − <sup>2</sup>)(*y* − <sup>5</sup>)<sup>2</sup>(*y* − <sup>7</sup>)<sup>3</sup>(*y* − <sup>11</sup>);

*Analytic*( *u*, *y*, *re* = −1..10, *im* = −2 ..10 );

7.00000000000000, 7.00000000000000, 7.00000000000000, 5.00000000000000, 5.00000000000000, 2.00000000000000, 11.00000000000000

According to Rouché's theorem, the denominator of Equation (6) has *b* zeros say, *γi* (*i* = 1, 2, ... , *b*) inside the unit circle. As *P*+*N* (*z*) converges in | *z* |≤ 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

$$P\_N^+(z) = T \left( \mathbb{C}(z) + \frac{f\_1(z)}{\prod\_{\substack{b+s\\i=b+1}}^{b+s} (z-\gamma\_i)} \right)' \tag{7}$$

where *<sup>C</sup>*(*z*) = *n*0 ∑ *i*=0 *ciz<sup>i</sup>* and *T* is a normalizing constant. Note that when 2*m* < *s*, then *<sup>C</sup>*(*z*) 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.


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

$$P\_N^+(z) = T\left(\mathbb{C}(z) + \sum\_{i=b+1}^{b+s} \frac{M\_i}{z - \gamma\_i}\right),\tag{8}$$

where

$$M\_i = \frac{f\_1(\gamma\_i)}{\prod\_{j=b+1, j \neq i}^{b+s} (\gamma\_i - \gamma\_j)}.$$

Using Equation (8), we have

$$P\_n^+ = \begin{cases} T\left(c\_n + \sum\_{i=b+1}^{b+s} \frac{-M\_i}{\gamma\_i^{n+1}}\right), & \text{if } n = 0, 1, \dots, n\_{0, \text{r}} \\\\ T\sum\_{i=b+1}^{b+s} \frac{-M\_i}{\gamma\_i^{n+1}}, & \text{if } n = n\_0 + 1, n\_0 + 2, \dots, N. \end{cases} \tag{9}$$

Employing the normalization condition *N* ∑ *<sup>n</sup>*=0 *<sup>P</sup>*+*n* = 1, we obtain the only unknown *T* as

$$T = \begin{cases} \left(\sum\_{i=0}^{n\_0} c\_i - \sum\_{i=b+1}^{b+s} \frac{M\_i}{\gamma\_i} \times \frac{1 - \gamma\_i^{-(N+1)}}{1 - \gamma\_i^{-1}}\right)^{-1}, & \text{if } \rho \neq 1, \\\\ \left(\sum\_{i=0}^{n\_0} c\_i - M\_{b+1}(N+1) - \sum\_{i=b+2}^{b+s} \frac{M\_i}{\Psi\_i} \times \frac{1 - \Psi\_i^{-(N+1)}}{1 - \Psi\_i^{-1}}\right)^{-1}, & \text{if } \rho = 1. \end{cases} \tag{10}$$

Thus, once all the roots are known, we can ge<sup>t</sup> the distribution for the number in queue.

**Remark 6.** *It may be noted that it is also possible to find the probabilities* {*P*+*n* }*N*0 *by assuming the solution of the form P*+*j* = *Czj, where C* = 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 postdeparture epochs by basic probabilistic reasoning and discrete-time renewal theory. In steady-state, let {*Pj*}*<sup>N</sup>*0 and {*<sup>P</sup>*<sup>−</sup>*j* }*N*0 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 *Pj* = *P*−*j* , ∀ *j* = 0, 1, ... , *N*; for details, see Boxma and Groenendijk [27]. If the server is idle, there are < *a* jobs in the queue. Suppose *Nq* is the number of tasks in the queue at some random time. At a random epoch, the steady-state probabilities are *Pn*,<sup>0</sup> = *P*(*Nq* = *n*, server idle), 0 ≤ *n* ≤ *a* − 1, and *Pn*,<sup>1</sup> = *P*(*Nq* = *n*, server busy), 0 ≤ *n* ≤ *N*. Given this,

$$P\_j = \begin{cases} P\_{j,0} + P\_{j,1} & \text{if } 0 \le j \le a - 1 \\ P\_{j,1} & \text{if } a \le j \le N \end{cases}$$

Let the limiting pmf of the elapsed service time be *s*ˆ-, which is determined by *s* ˆ - = *μ* ∞ ∑ *<sup>r</sup>*=-+1 *sr*, - ≥ 0 (see, [5] (p. 20)), and ˆ *k*- be the probability that the number of arrivals during an elapsed service time is -. This yields

$$\pounds\_{\ell} = \sum\_{i=\ell}^{\infty} \binom{i}{\ell} \lambda^{\ell} (1-\lambda)^{i-\ell} \pounds\_{i\prime} \; \ell = 0, 1, \dots$$

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,

$$E^\* = E(\mathcal{S}) \left( 1 - \sum\_{i=0}^{a-1} P\_i^+ \right) + \sum\_{i=0}^{a-1} P\_i^+ \left( \frac{a-i}{\lambda} + E(\mathcal{S}) \right) = E(\mathcal{S}) + \sum\_{i=0}^{a-1} P\_i^+ \frac{(a-i)}{\lambda} \cdot \mathcal{S}$$

**Remark 7.** *It can also be put down as E*∗ = *E*(*S*) + *a* ∑ *i*=1 *<sup>P</sup>*+*a*−*<sup>i</sup> iλ . When a* = 1*, it matches with Chaudhry and Chang [7].*

**Theorem 1.** *The random- and post-departure-epoch probabilities* {{*Pj*,<sup>0</sup>}*<sup>a</sup>*−<sup>1</sup> 0 , {*Pj*,<sup>1</sup>}*<sup>N</sup>*<sup>0</sup> }*, and* {*P*+*j* }*N*0 *are related by*

$$P\_{j,0} = \frac{\sum\_{i=0}^{j} P\_i^+}{\sum\_{i=0}^{a-1} (a-i)P\_i^+ + \lambda E(S)}, \quad 0 \le j \le a-1,\tag{11}$$

$$\begin{split} P\_{j,1} &= \left(1 - \sum\_{i=0}^{a-1} P\_{i,0}\right) \cdot \left[\sum\_{i=0}^{a-1} P\_i^+ \hat{k}\_j + \sum\_{i=a}^{b-1} P\_i^+ \left(\sum\_{m=i}^b y\_m\right) \hat{k}\_j + \sum\_{i=a+1}^{b-1} P\_i^+ \sum\_{m=a}^{i-1} y\_m \hat{k}\_{j-i+m} \\ &+ \sum\_{i=b}^N P\_i^+ \sum\_{m=a}^b y\_m \hat{k}\_{j-i+m} \right], \ 0 \leq j \leq N-1. \end{split} \tag{12}$$

*Finally, PN*,<sup>1</sup> *can be found from PN*,<sup>1</sup> = 1 − *<sup>a</sup>*−1 ∑ *j*=0 *Pj*,<sup>0</sup> − *N*−1 ∑ *j*=0 *Pj*,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 (*Pidle*). Let *E*(*I*) be the mean idle period. Using the definition of *E*(*S*) and *<sup>E</sup>*(*I*), we have

 $P\_{\text{idle}} = \frac{E(I)}{E(I) + E(S)},$  
$$P\_{j,0} = P\_{\text{idle}} \times P(\text{fraction of idle period}) = \frac{E(I)}{E(I) + E(S)} \times \frac{\frac{1}{\lambda} \sum\_{i=0}^{j} P\_i^+}{E(I)}, \ 0 \le j \le a - 1,$$

where *E*(*I*) = 1*λ <sup>a</sup>*−1 ∑ *i*=0 *P*+*i* (*a* − *i*). We employ system state conditioning and discreet renewal theory to find *Pj*,1. The processor is active with probability (1 − *<sup>a</sup>*−1 ∑ *i*=0*Pi*,<sup>0</sup>); thus,

$$P\_{j,1} = P(N\_q = j \text{,processor active}) = (1 - \sum\_{i=0}^{a-1} P\_{i,0})P(N\_q = j \mid \text{processor active}) \tag{13}$$

Assuming that *kej* is the number of jobs that come following an embedded Markov point until the elapsed service time, we have

$$P(N\_q = j \mid \text{processor active}) = \sum\_{i=0}^{a-1} P\_i^+ \sum\_{m=a}^b y\_m \hat{k}\_j + \sum\_{i=a}^{b-1} P\_i^+ \sum\_{m=i}^b y\_m \hat{k}\_j$$

$$+ \sum\_{i=a+1}^{b-1} P\_i^+ \sum\_{m=a}^{i-1} y\_m \hat{k}\_{j-i+m} + \sum\_{i=b}^N P\_i^+ \sum\_{m=a}^b y\_m \hat{k}\_{j-i+m},\tag{14}$$

$$0 \le j \le N - 1. \tag{14}$$

Putting together (13) and (14), we obtain (12). We can obtain *PN*,<sup>1</sup> using normalization condition. Thus, we obtain random epoch probabilities {{*Pj*,<sup>0</sup>}*<sup>a</sup>*−<sup>1</sup> 0 , {*Pj*,<sup>1</sup>}*<sup>N</sup>*<sup>0</sup> } in connection with post departure epoch probabilities {*P*+*j* }*N*0 .

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 queuelength distribution exactly before arrival of job will be equal to that of *Pj*,<sup>0</sup> and *Pj*,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 (*Lq*) is given by

$$L\_q = \sum\_{j=0}^{a-1} j \, \, P\_{\bar{j},0} + \sum\_{j=0}^{N} j \, \, P\_{\bar{j},1}$$

The probability of the processor being busy (PB) at some random moment is specified by 1 − ∑*<sup>a</sup>*−<sup>1</sup> *j*=0 *Pj*,0. Due to the BASTA property, the loss or blocking probability (*PBL*) is given by *PBL* = *PN*,1. Since the effective arrival rate *λe* = *λ*(1 − *PN*,<sup>1</sup>), we can obtain the average wait time in the queue (*Wq*) by employing Little's law as *Wq* = *Lq λe* . The reported result of *Wq* in [13] is incorrect. They have applied the effective arrival rate as *λ* instead of *λ<sup>e</sup>*.

**Remark 8.** *If ρ* < 1 *and N* → <sup>∞</sup>*, then λe* = *λ and Lq* = *<sup>a</sup>*−1 ∑ *j*=0 *j Pj*,<sup>0</sup> + ∞ ∑ *j*=0 *j Pj*,1*. Using Little's law, the average waiting time in the queue (Wq) can be computed as Wq* = *Lqλ .*

#### **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 (*Lq*), the average waiting time in the queue (*Wq*), and the probability of loss (*Ploss*). 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*/*NB<sup>a</sup>*,*Y*2 /1/10*queue. Consider the distribution of service time as being in two stage negative binomial distribution (NB) with PGF <sup>S</sup>*(*z*) = *μz* <sup>1</sup>−*μ*¯*z*<sup>2</sup>*. Here, we consider the same parameters as in Table 1 of the paper [13] to compare the results. The parameters are E*(*S*) = 5, *y*2 = 0.2, *y*3 = 0.1*, y*4 = 0.7*, and E*(*Y*) = 3.5*. The arrival rates are* 0.14, 0.35*, and* 0.56 *with corresponding traffic intensities ρ* = 0.2, 0.5*, and* 0.8*, respectively. To show the evaluation process, let us assume λ* = 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^+(z) = 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* (*z*) *has two roots outside the unit circle when arrival rates are* 0.35 *and* 0.56*, and we have from Equation* (8)*,*

$$P\_N^+(z) = T\left(c\_0 + \frac{M\_{b+1}}{z - \gamma\_{b+1}} + \frac{M\_{b+2}}{z - \gamma\_{b+2}}\right),$$

*where γb*+<sup>1</sup> = 2.106094*, γb*+<sup>2</sup> = 3.454495*, T* = 6.256038*, c*0 = 0.053212, *Mb*+<sup>1</sup> = <sup>−</sup>0.619857*, Mb*+<sup>2</sup> = 1.113392 *and γb*+<sup>1</sup> = 1.304299*, γb*+<sup>2</sup> = 2.688401*, T* = 6.867372*, c*0 = 0.025779*, Mb*+<sup>1</sup> = <sup>−</sup>0.081358*, Mb*+<sup>2</sup> = 0.224798*, respectively.*

Now, we can find post-departure epoch probabilities from Equation (9). The results are presented in Table 1. 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.


**Table 1.** Queue-length distributions at various epochs for the *Geo*/*NB*2,*Y*2 /1/10 queue.

**Remark 9.** *It may be noted that P*+*N* (*z*) *is a polynomial in both cases, as can be seen in Table 1. The same applies to other cases as well.*

**Example 2.** *Geo*/*DPHa*,*<sup>Y</sup>*/1/*Nqueue. The service-time distribution is assumed to be discrete phase-type (DPH) having sk* = *<sup>α</sup>Tk−***1***T***0***, k* = 1, 2, ...*, T***<sup>0</sup>** = *e* − *T<sup>e</sup>, where e is the appropriate column vector with all elements equal in size. This gives the PGF of service-time distribution as <sup>S</sup>*(*z*) = *zα***(***<sup>I</sup> − zT***)***−***1***T***0***,* |*z*| ≤ 1*. Table 2 shows the queue length distributions at different times employing the DPH service time distribution. For the first example of Table 2, we suppose* 0.200.05

*α* = \$0.40 0.50 0.10%*,* **T** = ⎡⎣0.10 0.30 0.15 0.10 0.20 0.50 0.10⎤⎦*,*

*with E*(*S*) = 2.005267 *and the other parameters are λ* = 0.7, *N* = 15, *y*3 = 0.5, *y*5 = 0.3, *y*8 = 0.2 *with E*(*Y*) = 4.6*. Here the denominator of P*+*N* (*z*) *has eleven distinct roots, out of which three roots are outside the unit circle, and they are γb*+<sup>1</sup> = 2.254374*, γb*+<sup>2</sup> = <sup>−</sup>9.788249*, and γb*+3 = <sup>−</sup>222.447161*. From Equation* (8)*,*

$$P\_N^+(z) = T \left( -3.63007 + \frac{1.301912}{z - 2.254374} - \frac{0.702111}{z + 9.788249} + \frac{930.333318}{z + 222.447161} \right), \tag{15}$$

*where T* = <sup>−</sup>1.755933*.*

*In the second example of Table 2, we assume α* = \$0.60 0.40%*, and T* = 0.5 1/3 1/3 1/3*with E*(*S*) = 4.2 *and the other parameters are λ* = 0.84, *N* = 50, *y*3 = 0.7, *y*4 = 0.2, *y*5 = 0.1 *with E*(*Y*) = 3.4*. As ρ* = 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)*,*

$$P\_N^+(z) = 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^{-19}}{z - 16.100784} \right), \tag{16}$$

*where T* = <sup>−</sup>7.569441*.*


**Table 2.** Queue-length distributions at various epochs for the *Geo*/*DPH*3,*<sup>Y</sup>*/1/15 and *Geo*/*DPH*3,*<sup>Y</sup>*/1/50 queues.

**Remark 10.** *It may be noted that Equations* (15) *and* (16) *are polynomials since the coefficients of zb 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*/*MGeo*2,*<sup>Y</sup>* 2 /1/10 *queue. Here, the service time is a mixture of two geometric distributions with PGF <sup>S</sup>*(*z*) = *ς*1 *μ*1*<sup>z</sup>* <sup>1</sup>−(<sup>1</sup>−*μ*1)*<sup>z</sup>* + *ς*2 *μ*2*<sup>z</sup>* <sup>1</sup>−(<sup>1</sup>−*μ*2)*<sup>z</sup> , where ς*1 + *ς*2 = 1*. The parameters are taken as E*(*S*) = 2.83333, *N* = 10, *ς*1 = 0.6, *μ*1 = 0.4, *μ*2 = 0.3, *y*2 = 0.2, *y*3 = 0.1, *y*4 = 0.7, *E*(*Y*) = 3.5*, λ* = 0.14*. The denominator of P*<sup>+</sup> *N* (*z*) *has six distinct roots, out of which two roots are outside the unit circle, and they are γb*+<sup>1</sup> = 4.031316 *and γb*+<sup>2</sup> = 5.727647*. From Equation* (8)*,*

$$P\_N^+(z) = T \left( -0.561775 - \frac{1.819848}{z - 4.031316} - \frac{4.545280}{z - 5.727647} \right) / 1$$

*where T* = 1.0000001*.*

*(ii) Geo*/*D*2,*<sup>Y</sup>*/1/50 *queue. Here, we consider service-time distribution as deterministic with PGF <sup>S</sup>*(*z*) = *zk for some constant k* = 4,*sk* = 1*. The parameters are taken as E*(*S*) = 4, *N* = 50, *y*2 = 0.5, *y*3 = 0.2, *y*4 = 0.1, *y*5 = 0.1, *y*6 = 0.1, *E*(*Y*) = 3.1*, and λ* = 0.4*. In this case, the denominator of P*<sup>+</sup> *N* (*z*) *has eight distinct roots, out of which two roots are outside the unit circle, and they are γb*+<sup>1</sup> = 5.019701 *and γb*+<sup>2</sup> = <sup>−</sup>11.795711*. From Equation* (8)*,*

$$P\_N^+(z) = T \left( -1.316734z^2 + 1.021799z - 102.664691 + \frac{879.864750}{z + 11.795711} - \frac{141.482007}{z - 5.019701} \right),$$

*where T* = 1*. Table 3 presents queue-length distributions at various epochs when the service-time distributions are a mixture of geometric and deterministic.*


**Table 3.** Queue-length distributions at various epochs for the *Geo*/*MGeo*2,*<sup>Y</sup>* 2 /1/10 and *Geo*/*D*2,*<sup>Y</sup>*/1/50 queues.

**Example 4.** *Here, we consider two cases. Table 4 presents the results of Geo*/*Ga*,*<sup>Y</sup>*/1/∞ *which can be found from the Geo*/*Ga*,*<sup>Y</sup>*/1/*<sup>N</sup> system by assuming ρ* < 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 ρ* < 1 *and N sufficiently large.*

*(i) Geo*/*NBa*,*<sup>Y</sup>*/1/∞ *queue. Here, we assume negative binomial (NB) service time distribution and the parameters are taken as λ* = 0.703, *E*(*S*) = 5, *y*2 = 0.5, *y*4 = 0.2, *y*5 = 0.1, *y*6 = 0.1, *y*8 = 0.1 *with E*(*Y*) = 3.7*. From Equation* (8)*, we have*

$$P\_N^+(z) = T \left( 0.005129 + \frac{0.03006}{z - 2.520637} - \frac{0.00745}{z - 1.042657} \right) \text{s}$$

*where T* = 6.250024*. (ii) Geo*/*DPHa*,*<sup>Y</sup>*/1/∞ *queue. For a DPH service time distribution, the settings are chosenas α* = \$0.60 0.40% *and* **T** = 0.5 1/3 1/3 1/3 *.*

*with E*(*S*) = 4.2 *and the other parameters are λ* = 0.58, *y*2 = 0.5, *y*3 = 0.3, *y*4 = 0.2*, and E*(*Y*) = 2.7*. From Equation* (8)*, we have*

$$P\_N^+(z) = 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),$$

*where T* = 20.250*.*


**Table 4.** Queue-length distributions at various epochs when *N* → ∞.

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 *ρ* = 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 2.** Time (in seconds) needed to calculate post-departure probabilities.

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 *λ* = 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 *ρ* = 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{aligned} z^{36} \big(-0.886 + 0.486z\big)^{4} - 0.0256(0.19 + 0.81z)^{4} \big(0.1 + 0.4z^{32} + 0.1z^{30} \\ + 0.1z^{24} + 0.2z^{11} + 0.1z^{6}\big). \end{aligned} \tag{17}$$

**Figure 3.** The 40 roots of Equation (17) when NB service-time distribution.

## **6. Conclusions**

This article focuses on the *Geo*/*Ga*,*<sup>Y</sup>*/1/*<sup>N</sup>* 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).

**Author Contributions:** Conceptualization, M.C. and V.G.; methodology, V.G. and M.C.; Writing— original draft preparation, V.G.; supervision, M.C.; funding acquisition, M.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by FMAS IO 757193 (C143).

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. The Continuous-Time Case**

Here, we consider the relation between the discrete-time *Geo*/*Ga*,*<sup>Y</sup>*/1/*<sup>N</sup>* queue and its continuous-time analogue. Let the time axis be divided into periods of uniform length Δ*u* with Δ*u* > 0 sufficiently small. In *Geo*/*Ga*,*<sup>Y</sup>*/1/*N*, since the inter-arrival times (*u*) follow geometric distribution, the arrivals will follow binomial distribution which, as we have seen earlier, leads to PGF for binomial distribution. In the continuous-time case, geometric tends to exponential, and binomial tends to Poisson distribution, and the PGF tends to

Laplace transform. Let us discuss this analytically. Assume that the inter-arrival times in the case of *<sup>M</sup>*/*Ga*,*<sup>Y</sup>*/1/*<sup>N</sup>* have a rate *α*. Then, *λ* = *α*Δ*u* + *<sup>o</sup>*(<sup>Δ</sup>*u*). In the discrete-case, let the service times *S* be in multiples of Δ*u* with probability *P*(*S* = -<sup>Δ</sup>*u*) = *s*- and ∞ ∑ -=1 *s*- = 1. Further, let *n*Δ*u* = *vn*, where the interval [0, *vn*] is divided into *n* intervals of length Δ*u*. The PGF of an arrival (or no arrival) in the interval (*<sup>v</sup>*-, *v*-+<sup>1</sup>) is (1 − *λ* + *<sup>λ</sup>z*). If we denote the probability density function of service times by *h*(·), then

$$\begin{aligned} P(\text{service finishes in } (\iota, \iota + \Delta \iota) \mid \text{service time} > \iota) &= h(\iota) \Delta \iota + o(\Delta \iota) \\ \text{and} \quad s\_\ell &= h(v\_\ell) \Delta \iota + o(\Delta \iota). \end{aligned}$$

When Δ → 0, the PGF *S*(1 − *λ* + *<sup>λ</sup>z*) changes to a Laplace transform. Using the definition of Lebesgue integration and taking the limit as Δ → 0 and *λ* = *<sup>α</sup>*Δ, we have

$$\lim\_{\Delta \to 0} K(z) = \lim\_{\Delta \to 0} \sum\_{\ell=1}^{\infty} s\_\ell (1 - \lambda + \lambda z)^\ell = \bar{h}(\alpha - \alpha z).$$

The proof of the above is not discussed in detail here since the method applied can be found in [28]. Now, using *λ* = *<sup>α</sup>*Δ, *<sup>K</sup>*(*z*) = ¯ *h*(*α* − *αz*) in (5) and taking the limit as Δ → 0, we have

$$P\_N^+(z) = \frac{\bar{h}(a - az) \left[ \sum\_{i=0}^{a-1} P\_i^+ \left( z^b - z^i Y'(z) \right) + \sum\_{i=a}^{b-1} P\_i^+ \left( \sum\_{r=i}^b y\_r z^b - z^i \sum\_{r=i}^b y\_r z^{b-r} \right) \right]}{z^b - \bar{h}(a - az) Y'(z)},$$

the connections for *<sup>M</sup>*/*Ga*,*<sup>Y</sup>*/1/*<sup>N</sup>* system. Taking *a* = 1, we have

$$P\_N^+(z) = \frac{h(\alpha - \alpha z) \left[ \sum\_{i=0}^{b-1} P\_i^+ z^b \left( \sum\_{r=i}^b y\_r - z^i \sum\_{r=i}^b y\_r z^{-r} \right) \right]}{z^b - \tilde{h}(\alpha - \alpha z) Y'(z)},$$

which matcheswith [29].

Note that we cannot obtain results of discrete-time analogue from [29], that is, the converse is not true.
