A Novel Computational Procedure for the Waiting-Time Distribution (In the Queue) for Bulk-Service Finite-Buffer Queues with Poisson Input

: In this paper, we discuss the waiting-time distribution for a ﬁnite-space, single-server queueing system, in which customers arrive singly following a Poisson process and the server operates under ( a , b ) -bulk service rule. The queueing system has a ﬁnite-buffer capacity ‘ N ’ excluding the batch in service. The service-time distribution of batches follows a general distribution, which is independent of the arrival process. We ﬁrst develop an alternative approach of obtaining the probability distribution for the queue length at a post-departure epoch of a batch and, subsequently, the probability distribution for the queue length at a random epoch using an embedded Markov chain, Markov renewal theory and the semi-Markov process. The waiting-time distribution (in the queue) of a random customer is derived using the functional relation between the probability generating function (pgf) for the queue-length distribution and the Laplace–Stieltjes transform (LST) of the queueing-time distribution for a random customer. Using LSTs, we discuss the derivation of the probability density function of a random customer’s waiting time and its numerical implementations.


Introduction
Single-server queueing models have been extensively studied in the literature.However, many practical situations exist where customers need to be served in batches.Queueing models of this type are called general bulk-service (GBS) queues.Examples of bulkservice queues can be seen in the unloading and loading of cargoes at a sea port, in traffic-signal systems, and in a mass transportation system.The GBS rule is stated as follows.If n (0 ≤ n < a) customers are waiting after the completion of a batch service, then the server has to wait until the number reaches 'a customers.However, the entire queue is served if the number of waiting customers in the queue is greater than or equal to a but less than or equal to b.The first 'b customers will be taken into service if the number present in the queue is greater than b.This rule has been widely investigated in the literature, e.g., Chaudhry and Templeton [1] and Dshalalow [2,3].An excellent survey of various types of infinite-buffer queues with bulk arrival, bulk service and vacations can be found in Dshalalow [4].In general, the analysis of queueing models of this type is difficult, since their probability generating functions involve unknown probabilities in the numerators.Neuts [5] discussed the bulk-service queue using the theory of the semi-Markov process and obtained queue-length distribution.Borthakur [6] studied the bulk-service queue and derived the stationary probability distribution for the queue length directly as the sum of a constant multiplied by geometric terms, where the common ratios are the root of the associated characteristic equation.Following this, Medhi [7] obtained the waiting-time density function for M/M (a,b) /1 queue in terms of roots.Later, Easton and Chaudhry [8] extended these results and discussed the waiting-time distribution for the E k /M (a,b) /1 queueing system.We can find an elaborate presentation of bulk-service queues in the books by Chaudhry and Templeton [1] and Medhi [9], where a complete list of references up to the early 1980s.One may note that it has been reported in the literature that the analysis of finite-buffer queueing models by solving simultaneous linear equations may pose some problems for some particular types of queueing models; see the remarks made by Powell ([10], pp.64, 141).Chaudhry et al. [11] discussed the computational analysis of the M/G (a,b) /1 and related non-bulk queues.Brière and Chaudhry [12] analyzed a batchservice queue with a random service batch size, i.e., M/G Y /1 queue.They obtained the probability distribution for the queue length at different epochs by finding the roots of the associated characteristic equation.Waiting-time analysis in a finite-buffer M/G/1 queueing model was first investigated by Finch [13] in his pioneering work in this area of research.However, no analytical results have been reported in the literature on queues to analyze the waiting-time distribution of a finite-buffer M/G/1-type bulk-service queue.A few results related to stationary queue-length distributions in a finite-buffer single-server queue using the roots' method are addressed in the literature, e.g., see, Chaudhry et al. [14] and Chaudhry and Goswami [15].Further, analysis of the finite-buffer bulk-service M/G Y /1 queue was performed by Singh [16] and queue-length distributions were obtained.Later, Chaudhry and Gupta [17] performed the analytic analysis of the finite-buffer bulk-service M/G (a,b) /1/N queueing model and obtained the queue-length distributions at a postdeparture and arbitrary epochs using the supplementary variable technique.However, no analytical results have been reported in the literature related to the analysis of the waiting-time distribution in a finite-buffer M/G/1-type bulk-service queue using roots.
We encounter finite-buffer bulk-service queueing systems in many real-world scenarios; see Gold and Tran-Gia [18] for applications related to manufacturing processes.In telecommunications networks, information units are stored in the system if a server is busy.If an arrival finds the buffer space full, then the arriving customer is blocked and considered to be lost.The main interest of the system designer is to provide sufficient buffer space so that the probability of blocking or loss is small.Thus, the blocking/loss probabilities are an important measure of concern for a system designer.In cloud computing, transmissions of user requests are of various types (e.g., data, voice, video, images).There are many applications which need bulk data to be transferred in the web which are accessed by a large number of clients at the same time.As per the service level agreement (SLA) all the requests received by the clients may be processed in bulk.The instances of the target web application running into virtual machines (VMs) act as service centers to process the requests in the queue.We figured out the process of client requests in a single web application with the following features: The inter-arrival times between two successive client requests may be identified by a Poisson process.When bulk requests are submitted in the queue, the requests are forwarded to the VM that is currently idle because of the lack of a sufficient number of requests.The requests are served in batches with a minimum of a and a maximum of b requests for this VM, where 1 ≤ a ≤ b ≤ N. The system under consideration contains multiple homogeneous VMs which render service in order of task request arrivals, i.e., first come, first served (FCFS).For a recent survey of the system model of bulk service on cloud, the readers are referred to Goswami et al. [19].In recent years, considering batch-size-dependent service processes, the finite-buffer capacity /1 queueing systems have been studied by Banerjee et al. [20], and Pradhan and Gupta [21], respectively.This paper carries out the analytic analysis of the finite buffer M/G (a,b) /1/N queue using the method of the Markov renewal theory, the embedded Markov chain, the semi-Markov process and the roots.As stated earlier, the service to the queueing system is provided in batches of minimum size a and maximum size b (1 ≤ a ≤ b).We obtain the steady-state, queue-length distribution at post-departure and random epochs.As stated above, using the functional relationship, we derive the probability density function for a random customer's waiting-time distribution (in the queue and in the system).Finally, to extend the work done by Chaudhry and Gupta [17] for the well-known model M/G (a,b) /1/N, which deals with the number in the queue, this paper first gives a simple procedure to get the probability generation function (pgf) of a number in the queue by using the method of the Markov renewal theory, the embedded Markov chain and the semi-Markov process.Using the derived pgf, we deal with the probability density function for the waiting-time distribution (in queue as well as in the system) of a random customer.In this single-server queueing model, the arrival process of the customers is assumed to be a Poisson process with an average rate, λ.The queued customers are served in batches of minimum size, 'a', with a maximum threshold size, 'b'; however, if there are fewer than 'a' customers, the server must wait until the queue accumulates 'a' customers.The number 'N' denotes the buffer capacity of the queueing system, excluding the size of a batch with the server, and it is assumed that 1 ≤ a ≤ b ≤ N. The service times of batches are independent identically distributed random variables (i.i.d.r.vs.) denoted by the notation S, where S has the probability density function (pdf) b(t), cumulative distribution function (DF) B(t), Laplace-Stieltjes transform (LST) B * (s) and mean service time

Brief
Let us define the conditional probability k j (j ≥ 0), which denotes the probability that j customers arrive in the queueing system during a service time of a batch.Then, The probability generating function (pgf) of {k j } ∞ 0 is denoted by K(z).It can be easily seen that K(z) = B * (λ − λz); see Keilson and Servi [22] for a detailed derivation.Similarly, we denote the conditional probability k j (j ≥ 0), which represents the probability that j customers arrive in the queueing system during a stationary elapsed service time of a batch.This leads to: Similar to the above result, if we define the pgf of { k j } ∞ 0 by K(z), then it can be seen that . The probability distribution for the queue length at a random epoch can be obtained through relations among post-departure and random epoch probabilities using the embedded Markov chain, the Markov renewal theory and the semi-Markov processes.For this, consider the queue at the batch service completion epochs.Let the successive batch service completion epochs be marked in the time-axis as t 0 , t 1 , t 2 , . . .The time epoch just after a service completion may be denoted by t + n .The state of the system at t + n is defined as {N q (t + n ), n ≥ 0}, where N q (t + n ) is the number of customers (0 ≤ N q (t + n ) ≤ N) in the queue after the n-th departure epoch of a batch.Thus, {N q (t + n ), n ≥ 0} becomes a Markov chain.This irreducible and aperiodic Markov chain becomes ergodic under any ρ >=< 1.In this discrete-time Markov chain, N q (t + n+1 ) is dependent on the random variable N q (t + n ) and another random variable A n+1 , which is independent of N q (t + n ).This gives: where x + = max{0, x} and A n+1 is the random number of arrivals during the (n + 1)th batch service time.When n → ∞, P(A n+1 = j) = k j (j ≥ 0).The various elements of the transition probability matrix (TPM) of this Markov chain are given by: where k c j = ∑ ∞ i=j k i (j ≥ 0).Thus, the limiting probabilities are defined as: Now p + (j) (0 ≤ j ≤ N) can be obtained by solving the system of equations P + = P + P, where , P is a (N + 1) × (N + 1) matrix whose elements are P ij , and after simplification of P + = P + P, we get: Here, it may be noted that we can find the pgf of {p + (n)} N 0 using (2) and (3).Multiplying (2) by z j and summing over j from 0 to N − 1 and adding (3) after multiplying with z N , we get after simplification the following expression: where p + * (z) = ∑ N n=0 p + (n)z n .Simplification of (4) by replacing ∑ N j=b+1 p + (j)z j = p + * (z) − ∑ b r=0 p + (r)z r in the right-hand side, we obtain: One can see that the numerator for the second term of Equation ( 5) has at least (N + b) degree in z.Since p + (n) (0 ≤ n ≤ N) is the coefficient of z n in the right-hand side of p + * (z) in Equation ( 5), then one may evaluate p + (n) using the following pgf: which is identical with the pgf of post-departure epoch probabilities in an infinite-buffer M/G (a,b) /1/∞ queueing model with the traffic intensity ρ < 1.The pgf of the infinitebuffer queueing model was mentioned earlier in the literature, e.g., see Equation (2.1) in [11] and also see Problem 11 in [1] (pp. 226-227)chaudhry1983first, where one may note that K(z) = B * (λ − λz).At this point, we intend to calculate the post-departure epoch probabilities in terms of the roots of the associated characteristic equation.Since rational LSTs (see Botta et al. [23]) cover a large number of distributions, we take K(z) = B * (λ − λz) = B * (s), where s = λ − λz as a rational function in 's'.If the LST of the service-time distribution has a non-rational LST, then the LST can be approximated by rational functions using the Padé approximation technique.
In view of this, we consider those service-time distributions that have rational LST of the form B * (s) = P(s)/Q(s), where the degree of the polynomial Q(s) is m and that of the polynomial P(s) is y (≤ m).Further, we make this assumption, as has been done by Botta et al. [23], since this assumption covers a large number of practical cases.Substituting in Equation ( 6), one can obtain the corresponding expression as given below, Since K(z) = f (z) g(z) , p + * (z) in Equation ( 7) becomes a rational function with the denominator (z b g(z) − f (z)) which is a polynomial of degree (b + m).This polynomial (z b g(z) − f (z)) has (b + m) zeros in the whole complex plane.This leads to the following characteristic equation (CE): which has (b + m) roots, namely, Remark 1.We could have avoided the derivation of pgf p + * (z), but since it gives a simpler procedure than what is available in the literature, we have included it here for completeness sake.This pgf is also used in deriving the waiting-time distribution.When software packages such as MAPLE and MATHEMATICA could not find (see Neuts [24], Kendall [25] and Kleinrock [26]) a large number of roots (they do now); Chaudhry responded to this (see Chaudhry [27,28]).In this connection, see also Chaudhry et al. [29].Here is just one example.It can be solved in many ways, but we are giving just one method using the package MAPLE: > restart : with(RootFinding) : 1., 2., 3., 3., 3., 3., 7., 7., 7.
The CE (8) has b + m roots.As the queueing model under consideration has a finite buffer, the following three cases arise depending on ρ: One may observe that as ρ (> 0) increases, from the remaining m roots of the CE (8), one positive real root gets closer to the origin from right to left.The above three cases have been discussed below.

•
We first discuss the case when ρ < 1.For this case, when the traffic intensity of the queueing model under consideration is less than one, i.e., ρ < 1, one can use Equation ( 6) to find the post-departure epoch probabilities as described below.The unknown probabilities {p + (r)} b−1 r=0 on the right-hand side of the pgf in Equation ( 7) can be obtained by using the roots of Equation ( 8), i.e., γ 1 = 1 and γ 2 , γ 3 , . . ., γ b .This leads to: and: Therefore, for ρ < 1, using the partial fraction p + * (z) can be written as follows: where D is the normalizing constant and l i 's are the constant coefficients of the partial fraction.Now, extracting the coefficients of z n from the right-hand side of ( 11), one can obtain the post-departure epoch probabilities as: where: One can obtain the unknown D in Equation ( 12) using the normalization condition as stated before and is given by: • For ρ = 1, we assume γ b+1 = 1, and in this case, Equations ( 11) and ( 12) are valid.The only exception is that l b+1 may be obtained as follows: and: It is also important to note that the normalization condition (10) does not work when the traffic intensity ρ = 1 as z = 1 is a repeated root of multiplicity two for the characteristic Equation ( 8).Thus, when ρ = 1, one may use Equations ( 16)- (18), which express the probabilities {p + (j)} N j=b in terms of the probabilities {p + (r)} b−1 r=0 and then one may use the normalization condition as p + * (1) = ∑ N j=0 p + (j) = 1.Alternatively, one can use the following normalization condition for ρ = 1: • Lastly, we discuss the case when ρ > 1.As p + * (z) is convergent/analytic inside |z| < γ b+1 , the first b roots of the characteristic Equation ( 8) inside the region |z| < γ b+1 should be the zeros of the numerator polynomial on the right-hand side of Equation ( 7).
The unknown probabilities {p + (r)} b−1 r=0 on the right-hand side of the pgf in Equation ( 7) can be obtained by the set of linear Equations (9).This set of homogeneous linear equations is not sufficient to give non-trivial solution for the unknown probabilities {p + (r)} b−1 r=0 and thus one must use normalization condition along with the above system of simultaneous Equations (9) in the following way.Thus, when ρ > 1, one should use Equations ( 2) and ( 3) by expressing the probabilities {p + (j)} N j=b in terms of the probabilities {p + (r)} b−1 r=0 as follows: and then one may use the normalization condition as follows: where the probabilities {p + (j)} N j=b must be expressed in terms of the unknown probabilities {p + (r)} b−1 r=0 using the above Equations ( 16)- (18).The solution of the above system of linear Equations ( 9) and (19) gives the unknown probabilities {p + (r)} b−1 r=0 and from the above expressions ( 16)-( 18) one can get the probabilities {p + (j)} N j=b in terms of the known probabilities {p + (r)} b−1 r=0 .It may be noted that the above procedure of finding post-departure epoch probabilities works when ρ <=> 1.One may remark here that the above proposed recursive scheme ( 16)-( 18) may sometimes be unstable, mainly when we are dealing with a very high buffer capacity (N → ∞) and due to the negative sign involved in the right-hand side of the recursive Equations ( 16) and (17).One may note here that it is possible to find post-departure epoch probabilities in closed form similar to the cases ρ <= 1 after getting the unknown probabilities {p + (r)} b−1 r=0 by using the above procedure.For this, we multiply both sides of Equation ( 7) by the factor (z − γ b+1 ) and obtain the following equation: Now, the left/right hand side of Equation ( 20) is convergent/analytic inside and on the unit circle |z| = 1.Thus, one may write Equation (20) in the following way: where q 0 , q 1 , and l i 's are the constant coefficients of the partial fraction.Equating the coefficient of the like powers of z in both sides of Equation ( 21) after a little algebraic calculation, we obtain: where: Let the coefficient of z 0 and z from the Taylors' series expansion of the right hand side of Equation ( 20) be c 0 and c 1 .Now, The last probability p + (N) can be found using the normalization condition (19) and is given by: p + (N) = 1 − ∑ N−1 j=0 p + (j).Thus, one can obtain post-departure epoch probabilities successfully for all possible values of offered load or traffic intensity ρ for this queueing model under consideration.Remark 2. It may be remarked here that all the probabilities p + (i) (0 ≤ i ≤ N) can also be computed by solving the system of linear equations p + P = p + (see [17], p. 98) using the GTH algorithm or the RG-factorization technique (extensions of the matrix-analytic methods, see Neuts [30] and Li [31]), where p + = [p + (0), p + (1), • • • , p + (N)] and P is the transition probability matrix of the Markov chain at the embedded post-departure epochs of a batch.In that case, the computational complexities of the GTH algorithm and the RG-factorization technique can be computed as O((N + 1) 4 ), (see [32]) and O(N 3 ) (see Li et al. [33]), respectively.On the other hand, in the case of the roots' method, which is adopted in this paper, the main goal is to calculate the m + b roots of the associated characteristic Equation ( 8) with the probability generating function given by Equation (7).Without a loss of generality, we assume that the degree of the characteristic Equation ( 8) is m + b.It should be mentioned here that we can numerically locate the roots of a polynomial using Newton's method.The root extraction using Newton's method and multiplication of two numbers has equal complexity, see [34].That is, if we want to calculate roots of the characteristic equation up to d 1 digits accuracy, then the computational complexity of per root extraction will be at most O(d 2  1 ).The next part is to evaluate b unknowns by solving b simultaneous linear equations.There are several methods used to solve simultaneous equations, e.g., Cramer's rule, Gaussian elimination, LU decomposition, etc.One of the most efficient methods to determine the b unknowns is LU decomposition with complexity O(b 3 ).The last part is a partial fraction decomposition which has computational complexity at most O(m 2 ) (see Xin [35]), where m is the degree of the denominator polynomial of the corresponding rational function.From our computational experience, we have seen that when the buffer capacity N becomes large and since b ≤ N, both the GTH algorithm and the RG-factorization technique take much longer to compute the post-departure epoch probabilities {p + (r)} N r=0 than the roots' method used in this paper.Further, the roots' method also gives explicit closed-form analytic results.In Section 5, we have presented a comparison of computation time for above discussed three methods.Lemma 1.The expression for the probability that the server is busy (ρ ) is given by: where T is mean inter-departure time of service batches.An expression for T may be obtained as follows: Proof.The server's state can be either busy or idle at any instant.As T is the mean interdeparture time of service batches, 1/T corresponds to the departure rate.Since T is the mean inter-departure time, it can be derived by considering two cases, (1) when the server is busy and (2) when the server is idle.Thus: which after simplification gives (26).The probability that the server is busy is given by Equation ( 25), which is the proportion of time the server is busy divided by the mean inter-departure time.

Probability Distribution for Queue Length at a Random Epoch
We develop relation between a distribution and the probability that the number of customers in the queue at the embedded Markov points (post-departure epochs) using the argument of the Markov renewal theory and the semi-Markov process.In terms of notations, the random or arbitrary epoch probabilities are defined as follows: where N q (t) is the number of customers in the queue at time t and ξ(t) is the server's state at time t.If at time epoch t the server is busy, then ξ(t) = 1, and if at time epoch t the server is idle with i ((0 ≤ i < a) customers waiting for service, then ξ(t) = 0. Therefore, p 1 (n) (0 ≤ n ≤ N) denotes the probability of n customers in the queue (excluding the size of a batch in service) when the server is busy and p 0 (n) (0 ≤ n ≤ a − 1) denotes the probability of n customers in the queue when the server is idle, where '1' stands for busy state of the server; and '0' stands for idle state of the server.
Proof.Using rate-in, rate-out, we can consider the state {N q (t) = n, ξ(t) = 0} in the limiting case of t → ∞.This can change by arrivals or service completions.Thus, we have: The left-hand side of Equation ( 33) is the rate of leaving out of the state {n, 0} by an arrival, while the first term of the right-hand side is the entering rate to the state {n, 0} after an arrival and the second term is the entering rate to the state {n, 0} after a departure.After a little simplification of (33), we get (30).Using rate-in, rate-out, we now consider the state {N q (t) = n, ξ(t) = 1} in the limiting case when t → ∞.This can change by service completions.As k j is the probability of number of customers that arrive following an embedded Markov point until stationary elapsed service time, then we have the following result: On the left-hand side of (34), we are leaving the state {j, 1} and on the right-hand side of (34), we are entering into the corresponding state {j, 1}.After simplifying (34), we obtain (31).An exactly similar argument as used in Equation ( 31) derives Equation (32).
One may note that the results derived in Equations ( 30) and ( 31) are actually the same as those presented by Chaudhry and Gupta [17]; see Appendix A for this derivation.The rate-in and rate-out principle used to derive Equations ( 30) and (31) are much simpler than the supplementary variable method used in [17].Moreover, one can check the correctness of the numerical results obtained by both the methods.Hereafter, let us define the pgf of {p 1 (j)} N 0 by P qb (z) = ∑ N j=0 p 1 (j)z j , which is the pgf of the random epoch probabilities {p 1 (j)} N 0 when the server is busy.Using (31) and (32), one can check that P qb (1) = ∑ N j=0 p 1 (j) = 1 µT , which is equal to ρ as derived in Lemma 2. One can check the fact that P qb (1) = 1 − ∑ a−1 i=0 p 0 (i) using Equations ( 30) and (25) in the following way: which is, in fact, true as 1 µT represents ρ which denotes the probability that the server is busy, i.e.: Moreover, it is interesting to note from Equation ( 35) that: Finally, one may note that using the fact that Poisson arrivals see time average (PASTA), the pre-arrival epoch probabilities are equal to the random epoch probabilities {p 0 (n) (0 ≤ n ≤ a − 1) and p 1 (j) (0 ≤ j ≤ N)}.Therefore, the probability of loss or blocking of a random arrival customer is given by P loss = p 1 (N).

Waiting-Time (In Queue) Distribution for a Random Customer
Let us denote as the waiting time in the queue of a random customer who belongs to those customers served in the n-th batch.Further, let W q (t) = lim n→∞ P(V in the limiting case as n → ∞, we denote V (n) q by V q .Let us define W * q (s) = ∞ 0 e −st dW q (t) with (s) ≥ 0. Now, following Chaudhry and Templeton [36], one can observe that during the busy or idle state of the server, the total probability of arrivals during the waiting time (in queue) of a random customer is equal to the total probability of possible number of customers waiting in the queue.Thus, considering the random nature of arrivals and busy/idle state of the server, we obtain: where λ (= λ(1 − P loss )) which is the customer's actual arrival rate due to finite-buffer capacity.We note that differentiating (37) with respect to z and setting z = 1, we obtain Little's law as follows: From Equation (37) using s = λ − λ z, after a little algebraic manipulation, we obtain: Now, using Padé approximation of the right-hand side of Equation ( 39), one may approximate W * q (s) by the rational function P(s) Q(s) . We assume the degree of the numerator P(s) as r 0 and that of the denominator Q(s) as r d (> r 0 ).One can check the correctness of Padé approximation using the fact that . Further, one can improve the Padé approximation using those two criteria.Since W * q (s) is convergent or analytic in (s) ≥ 0, the r d zeros of the denominator Q(s) may be assumed as s 1 , s 2 , . . ., s r d with (s i ) < 0. Using partial fractions, the right-hand side of Equation ( 39) may be written in the following form: Now, Equation ( 40) leads to the following equation: where , where j = 1, 2, . . ., r d .

Remark 3.
Here, we demonstrate how to practically approximate W * q (s) by the rational function using the Padé approximation through a numerical example.For this purpose, we consider PH-type service-time distribution with the representation (β β β, T T T), where β β β = (0.2, 0.8) and T T T = −2 0 0 −4 , so that µ = 3.333333.The other parameters are taken as a = 2, b = 5, λ = 10 with ρ = λ b•µ = 0.6.After computing p 0 (i)'s and p 1 (j)'s, the effective arrival rate of customers can be obtained as λ = λ(1 − P loss ) = 8.928835.From Equation (39) We have used the 'pade' command in MAPLE 2015 and approximated W * q (s) as P(s) where: Moreover, the correctness of the approximation can be checked by using the criteria given after the Equation (39): .000000, and − d ds All of the roots of Q(s) have real parts less than zero, which makes W * q (s) convergent in (s) > 0. Then we proceed further after writing in terms of partial fractions.Remark 4. One may note here that inversion of the LST presented in Equation ( 39) may be done using several other methods of inverting the LSTs (or transforms) reported in the literature (see the excellent reviews on the numerical transform inversion by Abate and Whitt [37] and Abate et al. [38]).Abate and Whitt [37] use the Fourier-series method to invert the generating functions and the LSTs numerically.The Fourier-series method involves numerically integrating a standard inversion integral employing the trapezoidal rule.The greatest difficulty in this case is approximately calculating the infinite series obtained from the inversion integral.However, the method described above is free of such an approximation of numerical integration.Moments: Moments of order k (≥ 1) for V q can be obtained using (41) as follows: k j e s j t dt, where (s where (s j ) < 0. Alternatively, one can get the k-th moment of the queueing-time distribution by differentiating (39) k times with respect to s and then setting s = 0 in the following way: which may work as a check on correctness of the moments obtained via (42).
Sojourn-time distribution: Sojourn-time distribution of a random customer is its queueing time, plus its batch service time.Therefore, the LST of sojourn time of random customer W * (s) is given by: Using (41), we can invert W * (s).It gives the pdf w(t) of the sojourn time distribution of a random customer and is given by: where (s j ) < 0. Now, using the probability density function w(t), one can get a moment of order k for the random variable V (sojourn time in the system for a random customer) as )dt which may be described as follows: The cumulative distribution function (CDF) of the queueing and sojourn times may be obtained by the following results: k j e s j u dt, t > 0, where (s j ) < 0,

Applications of the Queueing Model in the Performance Evaluation of Blockchain Systems
In a blockchain system, client transactions must wait for servers (e.g., miner, validator or orderer) to provide service (e.g., mining, validating or ordering) and finally get confirmed/served.One can find a detailed literature survey and a description of the performance evaluation of blockchain systems by Fan et al. [39].Different consensus processes of blockchain systems can be modelled as different types of queueing systems.In a queueing system, the most important part is the quantitative evaluation of its performance measures, such as what is the expected number of transactions in the system, what is the transaction throughput of the system, and what is the average waiting time in the queue or in the system (i.e., sojourn time).In the blockchain called Bitcoin [40], the ledger is maintained and updated by the mining process.In the mining process, a bunch of nodes, called miners, compete to solve challenging puzzle-like problems, which consume much computation power.Transactions issued by users are grouped into a container called a block.The mining competition winner who first finds the algorithmic puzzle answer specialized for the block has the right to add the new block to the blockchain and gets incentives accordingly.Kawase and Kasahara [41] first built a modified M/G B /1 queue (see Chaudhry and templeton [36] for the performance analysis of the M/G B /1 queueing model) with batch service to model the Bitcoin mining process.Li et al. [42] introduced a new blockchain queueing model by decomposing the mining process into two exponential service stages: block generation and blockchain-building processes.The sum of both stages' times is regarded as the transaction confirmation time, which a general service-time distribution may represent.The distribution of the key performance indicators, such as the the transaction waiting time, may be defined by the random amount of time a transaction has to wait before it gets confirmation; see Geissler et al. [43] for details.Geyer et al. [44] modelled the solo ordering process of Hyperledger Fabric as an M/M B /1 queueing system.To solve this model, the authors borrowed the results from the waiting time analysis of a general bulk service queueing model M/M (a,b) /1 (see Medhi [7]).They considered the batch size range as either a = b = B or a = 1 and b = B. Following the above discussions, we propose the queue length and waiting time analysis of a more general bulk-service M/G (a,b) /1/N queue to analyse the performance of a blockchain queueing system which includes the possibility that the buffer capacity N → ∞.See the following Figure 1.

Numerical Results
In this section, the validity of the results obtained in the previous sections has been tested through numerical experiments.Several numerical experiments have been carried out for different model parameters in this connection.Out of them, a few numerical examples are presented in this section.The experiments have been carried out using MAPLE 2015 for Windows 10 (64-bit) OS environment.The results have been presented in this paper up to six decimal places due to a lack of space.However, upon specific requests, the authors can provide high-precision numerical calculations exhibited through MAPLE coding.
Example 1.Here we will examine the case of traffic intensity (ρ) greater than 1.We construct a PH type distribution as the service-time distribution with the representation (β 3 β 3 β 3 , T 3 T 3 T 3 ) to serve that purpose, where β 3 The other parameters are chosen as λ = 20, N = 20, a = 3, b = 10 and ρ = 1.234667.Proceeding in a similar fashion mentioned in the above sections, we get the denominator of Equation ( 7) as: In this case, Equation ( 8) has eleven roots inside and on the unit disk and three roots outside the unit disk;  1, from which the loss probability of a random arrival customer can be obtained as P loss = p 1 (20) = 0.313256.
The obtained results are given in Table 5.From the table, the loss probability of a random arrival customer can be observed as P loss = p 1 (20) = 0.221596.The CDF of queueing-time distribution W q (t) for the model is obtained using the approach discussed in the paper as follows: W q (t) =0.999998 − 1.92626e −2.12386t cos(0.744301t)− 7.89988e −2.12386t sin(0.744301t) Moreover, using similar methodology conferred in the paper, the CDF of sojourn-time distribution W(t) for the model is obtained as follows:  Putting k = 1 and 2 in Equation ( 42), we obtain the following: Var(V q ) = 0.588934.
Similarly, a moment of order n can be obtained for the sojourn time in the system for a random customer by substituting k = 1 and k = 2 in Equation ( 46): Var(V) = 2.513934.
We have plotted loss probability (P loss ) and mean sojourn time (E[V]) against the buffer size N by varying it from 10 to 45 in Figure 5, from which it is apparent that as the system capacity rises, P loss decreases and the mean sojourn time increases linearly apparently.Example 6. Medhi [7] deals with the waiting-time distribution for the M/M (a,b) /1 queueing model under steady-state.In this example, we have analysed this model and compared it with the results of [7].Here, the service time is assumed to be exponentially distributed with a mean of 1 µ and B * (s) = µ µ+s .Let us take the parameters as a = 3, b = 6, N = 200, λ = 6.7, and µ = 1.7 with ρ = λ bµ = 0.656863.Following a similar procedure, as described in the examples mentioned earlier, we have obtained the following results and compared our results with that of [7] in Table 6.Moreover, we have calculated the CDF of the waiting time in the queue and compared it with the results of [7].Using the procedure discussed in the paper, W q (t) is obtained as follows: W q (t) =0.999999 − 0.0177191e −5.20699t − 0.125600e −5.10225t cos(10.5612t)− 0.257258e −5.10225t sin(10.5612t)− 0.856680e −0.905312t + 0.0000107034e −0.335552t , t > 0.
The results have been exhibited in Table 7 and Figure 6.Table 6.Performance measures of different quantities with that of [7].
In this example, we have analysed this model and compared it with [7]'s results.Here, the service time is assumed to be exponentially distributed with a mean of 1 µ = 1.0 and B * (s) = 1.0 1.0+s .Following a similar procedure, as described in the aforementioned examples, we have obtained the following results taking N = 300 and compared our results with that of [49] in Table 8.Also, we have compared them with the results of [7].The results have been exhibited in Table 8.Further, it can be seen that E(V) = 3.147639 exactly matches with that of [7] while E(V) = 3.147645 is given in [49].

Computation Time Comparison
Finally, the computation time for determining the queue-length distribution at a post-departure epoch of a batch in an M/PH (3,10) /1/N queueing system using the roots' method, the GTH algorithm and the RG-factorization technique (extensions of the matrixanalytic method) have been compared.The computation time for these three techniques (which are already discussed in the previous remark) has been presented graphically in Figure 7.We used MAPLE 2015 in a PC with Intel(R) Core i7 processor @3.40 GHz with 8 GB RAM and the Windows 10 environment for this experiment.The PH representation is taken the same as used for Example 1.The mean arrival rate (λ) is taken as λ = 20.The buffer capacity (N) varies from 11 to 25, and for each case, the computation time (in seconds) has been noted down for computing the stationary queue-length distribution at a post-departure epoch of a batch.The computation times against each buffer capacity are presented in Figure 7.Moreover, from Figure 7, it may be observed that the computation time using the roots' method is lower than that of the RG-factorization technique.Further, from Figure 7, one can also see that the computation time for the roots' method and GTH algorithm are almost the same in the given range of the buffer capacity.One may observe that the computation time using the roots' method and GTH algorithm is almost parallel to the x-axis in the given range of N, see Figure 7.

Conclusions
In this paper, the finite-buffer, bulk-service queue has been analysed, and the steadystate probability distributions at random-and post-departure epochs have been obtained.For this purpose, the roots of the characteristic equation are used.After that, the functional relation between the pgf of the queue-length distribution and the LST of the queueing-time distribution have been used to obtain the queueing time distribution when the server is busy.Moreover, we obtain the probability density functions of the sojourn time and the queueingtime distributions.The analysis of the stationary probabilities for the batch arrivals and batch service queues have problems in dealing with waiting-time distributions.Instead of the Poisson arrival process, one may consider a non-renewal (batch) arrival process (MAP or BMAP arrival processes) to the same bulk-service queue studied in this paper, and in such a case, the stationary distributions may be achieved using a solution procedure similar to the present article.The analysis of stationary probability vectors for the corresponding batch size dependent bulk-service service queues, such as BMAP/MSP

Figure 2 .
Figure 2. Effect of N on loss probabilities.

Figure 3 .
Figure 3.Effect of N on mean sojourn time.

Figure 4 .Example 5 .
Figure 4. Comparison among the CDF of the sojourn times for different values of N from Example 4 of Yu and Tang [49].Example 5.In this example, we consider phase-type (PH) service-time distribution with the

Figure 5 .
Figure 5. Buffer size (N) versus mean sojourn time (E[V]) and N versus loss probability (P loss ).

Figure 7 .
Figure 7. Computation time for different techniques with varying number of servers.
N (∞) and BMAP/G (a,b) n /1/N (∞) are interesting problems and may attract researchers' attention in the near future.

Table 7 .
Comparison between numerical results for CDF of queueing time of a random customer from Example 1 with the parameters λ = 6.7, µ = 1.7, a = 3, b = 6, N = 200 and ρ = 0.656863.Comparison of CDF of queueing time obtained from this paper against that ofMedhi [7].