**1. Introduction**

Queueing theory consists of a powerful tool for modelling and analytically studying many complex systems, such as computer networks, banks, telecommunications, manufacturing, and transportation systems. Compared to well-developed single-server non-bulk queueing systems, bulk-service systems have an extensive mathematical theory. They are more complex and harder to deal with. In a bulk-service queue, a group (or batch) of customers can be served simultaneously. Examples of their applications can be seen in shuttle-bus services, freight trains, express elevators, tour operators, and batch servicing in manufacturing processes. This topic, due to its perceived applicability, has attracted the attention of many researchers over several decades. At an early stage, some simple bulk-service models, such as single-server systems GI/M*<sup>b</sup>*/1 and M/M*a*/1 were studied by Shyu [1] and Gross et al. [2], respectively. Neuts [3] first introduced a quorum bulk service rule to create more complex models necessary to describe certain realistic situations. He considered a queueing system with Poisson arrivals and a general service-time distribution M/G*<sup>a</sup>*, *<sup>b</sup>*/1, where *a* is the quorum and *b* is the capacity of the server. Easton and Chaudhry [4] extended these results to the case where the inter-arrival times were Erlangian with the *η*-stage, <sup>E</sup>*η*/M*<sup>a</sup>*, *<sup>b</sup>*/1. Later, Chaudhry and Madill [5] gave a solution for a more general queueing system GI/M*<sup>a</sup>*, *<sup>b</sup>*/1. An alternate method was given in Neuts' book [6], wherein he describes the application of his matrix geometric approach to the GI/PH*<sup>a</sup>*, *<sup>b</sup>*/1 system, which has a phase-type service-time distribution. However, these systems are single-server queues. For many other variations of bulk-service queues, such

**Citation:** Chaudhry, M.; Gai, J. Analytic and Computational Analysis of GI/M*<sup>a</sup>*, *<sup>b</sup>*/*c* Queueing System. *Mathematics* **2022**, *10*, 3445. https://doi.org/10.3390/ math10193445

Academic Editors: Gurami Tsitsiashvili and Alexander Bochkov

Received: 16 August 2022 Accepted: 8 September 2022 Published: 22 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/).

as bulk service queues with vacations or bulk-service queues of the type M/G/1, one may view the survey paper written by Sasikala and Indhira [7]. In this survey, which had over 100 publications, most of the models considered were single server queues.

Multi-server queueing systems are an important class of queueing processes and have broad practical applications. However, such systems are more complex and harder to deal with compared to single-server queueing systems, especially when the inter-arrival time distribution is arbitrary. Medhi [8] investigated a queue with Poisson arrivals M/G*<sup>a</sup>*, *<sup>b</sup>*/*<sup>c</sup>*, but his method was not analytically tractable for *c* > 2. Related work has also been conducted by Sim [9] on M/M*<sup>a</sup>*, *<sup>b</sup>*/*c* by using algorithmic methods but no numerical results were given. Sim [10] solved the *η*-phase Erlangian arrivals <sup>E</sup>*η*/M*<sup>a</sup>*, *<sup>b</sup>*/*c* system for the random epoch probabilities in the steady state and discussed his results in the context of a transportation system. Adan and Resing [11] derived and presented the numerical results of the queue-length distributions for models M/COXIAN-2*<sup>a</sup>*, *<sup>b</sup>*/*c* and M/E*η a*, *<sup>b</sup>*/*<sup>c</sup>*. Compared to our model GI/M*<sup>a</sup>*, *<sup>b</sup>*/*<sup>c</sup>*, the most relevant model studied by other researchers was GI/M*<sup>b</sup>*/*<sup>c</sup>*, where the quorum was set to 1. Goswami et al. [12] solved the finite-buffer GI/M*<sup>b</sup>*/*c* model by the supplementary variable technique. Shyu [13], as well as Chaudhry and Templeton [14], dealt with the distribution of the number of customers in the system without considering the server being busy or idle. Therefore, there is no information regarding server utilization. Moreover, the numerical results for the system GI/M*<sup>b</sup>*/*c* are not available.

To make the model useful for applications, in this paper, we considered analytic and computational aspects to determine the performance of a complex bulk-service, multiserver queueing system GI/M*<sup>a</sup>*, *<sup>b</sup>*/*<sup>c</sup>*. The model GI/M*<sup>a</sup>*, *<sup>b</sup>*/*c* is an extension of the system GI/M*<sup>b</sup>*/*c* (Shyu [13] as well as by Chaudhry and Templeton [14]), by introducing quorum in the multi-server system GI/M*<sup>b</sup>*/*<sup>c</sup>*. A quorum refers to the minimum number of customers that are required in the waiting line before service commences, e.g., a ferry will not start until the quorum is met, or if we are dealing with transportation problems, a bus may not start until we have the quorum. This is an important policy desired by the service providers to reduce the business cost and maximize server utilization. The adding of the quorum policy makes the model closer to the real situation, but it also makes the model more complex to study. In view of this, a two-dimensional Markov chain has to be involved where the first dimension corresponds to the state of the servers (busy or idle) and the second dimension corresponds to the number of customers in the queue. We give an elegant analytic closed-form solution to obtain the queue-length distributions at three different epochs, such as pre-arrival epoch (p.a.e.), random epoch (r.e.), and postdeparture epoch (p.d.e.), not only for the system in a busy state, but also in an idle state. In the case of the idle state, the probabilities were obtained by simultaneously solving the *c* × *a* equations, some of which contained infinite series, which needed to be truncated to obtain the results. Instead of truncation, which leads to approximate results, we derived a closed-form solution and proposed an efficient algorithm to fix this problem. The model GI/M*<sup>a</sup>*, *<sup>b</sup>*/*c* that we considered includes most models ([1,2,4–6,8–10,13,14]) as special cases. Our model was validated in giving numerical results with the desired degree of accuracy and trivial computational costs. By selecting particular numbers for the parameters *a*, *b* and *c*, and inter-arrival time distributions, the numerical results produced by our model match the ones provided in those simpler models as expected.

The paper is organized as follows. In the following section, we describe the queueing model GI/M*<sup>a</sup>*, *<sup>b</sup>*/*<sup>c</sup>*, and establish a transition probability matrix (t.p.m.) for the system in Section 3. In Sections 4–6, we obtain the queue-length distributions at three different epochs, such as pre-arrival epoch (p.a.e.), random epoch (r.e.), and post-departure epoch (p.d.e.). To make the model useful for applications, sample numerical results are provided in Section 7.

#### **2. Model Description**

In this continuous-time queueing system GI/M*<sup>a</sup>*, *<sup>b</sup>*/*<sup>c</sup>*, there are *c* independent servers, each serving at the rate *μ*. The customers arrive at the rate *λ* according to a renewal process with an arbitrary inter-arrival time distribution *<sup>A</sup>*(*t*). One of the idle *c* servers starts the service as soon as the number of customers (including the new arriving customer) in the queue reaches quorum *a*. Each *c* server is able to serve up to *b* customers simultaneously. This indicates that if the server completes a service and finds less than the quorum *a* in the queue, it will become idle until *a* is reached. The service times of each server are independently–identically exponentially distributed random variables (i.i.e.d.r.v.s). We consider the system to be in a steady state with the traffic intensity *ρ* = *<sup>λ</sup>*/(*bcμ*) < 1. The queue discipline is first-come first-serve (FCFS) by batches.

#### **3. Transition Probability Matrix (t.p.m.)**

In the queueing system GI/M*<sup>a</sup>*, *<sup>b</sup>*/*<sup>c</sup>*, the states occurring at the instants immediately before the arrivals form an embedded Markov chain (I.M.C.). The state seen by an arriving customer can be described by (*Sn*, *<sup>n</sup>*), where *n* ≥ 0 is the queue-length and *Sn* is a supplementary flag defined as

$$S\_n = \begin{cases} I(k), & \text{if } k \text{ servers are idle}, \quad 1 \le k \le \mathfrak{c}, \quad 0 \le \mathfrak{n} \le \mathfrak{a} - 1, \\ B, & \text{if all servers are busy}, \quad \mathfrak{n} \ge 0. \end{cases}$$

We define the system as busy if all the servers are busy (*Sn* = *B*), and idle if at least one server is idle (*Sn* = *<sup>I</sup>*(*k*), *k* is the number of idle servers). The queue-length *n* can be written as *n* = *qb* + *n*0, 0 ≤ *n*0 ≤ *b* − 1, where *q* is the nearest lower non-negative integer of the fraction *<sup>n</sup>*/*b*, denoting the available number of full size batches (the batch size is *b*) in the queue waiting for service.

To build a t.p.m. of the system, we first define the following probabilities.

1. [*l*|*m*; *t*] and [*l*|*m*], where 0 ≤ *l* ≤ *m* ≤ *c*, and there are less than *a* customers waiting in the queue at the beginning of the period, thus *q* = 0. Here,

$$[l|m;t] = \binom{m}{l} (1 - e^{-\mu t})^l (e^{-\mu t})^{m-l}$$

is the conditional probability that *l* of *m* servers complete services during an interarrival period of duration *t*, given that *m* servers are busy (*c* − *m* servers are idle) at the beginning of the period. Moreover, [*l*|*m*] is defined as

$$\left[l\vert m\right] = \int\_0^\infty \left[l\vert m; t\right] dA(t), \qquad 0 \le l \le m \le c. \tag{1}$$

2. {*l*|*c*; *q*} is the conditional probability that *l* of *c* servers become idle during an interarrival period, given that all *c* servers are busy at the beginning of the period, and *q* (*q* ≥ 1) batches of customers are waiting for the services. Assume that a time *V* has elapsed when the last batch of *q* batches enters service. In this case, the *c* servers have been processed at a rate of *cμ* until time *V* has elapsed. When all *c* servers are busy, the number of departed batches follows a Poisson process with a rate *<sup>c</sup>μ*. The time *V* is Erlang-distributed, so it is the sum of *q* exponential random variables with a rate *<sup>c</sup>μ*, implying that the probability density function (p.d.f.) of *V* is given by

$$p(v) = \frac{(c\mu)(c\mu v)^{q-1}e^{-c\mu v}}{(q-1)!}, \qquad v > 0.$$

After all the waiting *q* batches leave the queue, there is time *t* − *V* remaining to have *l* batches processed. The probability that these *l* batches complete the service during period *t* − *V* is [*l*|*c*; *t* − *<sup>V</sup>*]. Therefore

$$\{I|c;q\} = \int\_0^\infty \int\_0^t \binom{c}{l} (1-e^{-(t-v)\mu})^l (e^{-(t-v)\mu})^{c-l} \frac{(c\mu)(c\mu v)^{q-1}e^{-c\mu v}}{(q-1)!} dv dA(t). \tag{2}$$

3. (*l*|*c*) is the conditional probability that *l* batches complete service during an interarrival period of duration *t*, given that all the *c* servers are busy at the beginning of the period and still busy at the end of the period. In this case, the number of batches served in time *t* is distributed as a Poisson process at a rate of *cμ*:

$$\langle l|c \rangle = \int\_0^\infty \frac{e^{-c\mu t} (c\mu t)^l}{l!} dA(t), \qquad l \ge 0. \tag{3}$$
