Analytic and Computational Analysis of GI/M^a,b/c Queueing System

Abstract

Bulk-service queueing systems have been widely applied in many areas in real life. While single-server queueing systems work in some cases, multi-servers can efficiently handle most complex applications. Bulk-service, multi-server queueing systems (compared to well-developed single-server queueing systems) are more complex and harder to deal with, especially when the inter-arrival time distributions are arbitrary. This paper deals with analytic and computational analyses of queue-length distributions for a complex bulk-service, multi-server queueing system GI/M^$a,b$/c, wherein inter-arrival times follow an arbitrary distribution, a is the quorum, and b is the capacity of each server; service times follow exponential distributions. The introduction of quorum a further increases the complexity of the model. In view of this, a two-dimensional Markov chain has to be involved. Currently, it appears that this system has not been addressed so far. An elegant analytic closed-form solution and an efficient algorithm to obtain the queue-length distributions at three different epochs, i.e., pre-arrival epoch (p.a.e.), random epoch (r.e.), and post-departure epoch (p.d.e.) are presented, when the servers are in busy and idle states, respectively.

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^b/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^$a,b$/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 $\eta$-stage, E_$\eta$/M^$a,b$/1. Later, Chaudhry and Madill [5] gave a solution for a more general queueing system GI/M^$a,b$/1. An alternate method was given in Neuts’ book [6], wherein he describes the application of his matrix geometric approach to the GI/PH^$a,b$/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 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^$a,b$/c, but his method was not analytically tractable for $c>2$. Related work has also been conducted by Sim [9] on M/M^$a,b$/c by using algorithmic methods but no numerical results were given. Sim [10] solved the $\eta$-phase Erlangian arrivals E_$\eta$/M^$a,b$/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^$a,b$/c and M/E_$\eta$^$a,b$/c. Compared to our model GI/M^$a,b$/c, the most relevant model studied by other researchers was GI/M^b/c, where the quorum was set to 1. Goswami et al. [12] solved the finite-buffer GI/M^b/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^b/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, multi-server queueing system GI/M^$a,b$/c. The model GI/M^$a,b$/c is an extension of the system GI/M^b/c (Shyu [13] as well as by Chaudhry and Templeton [14]), by introducing quorum in the multi-server system GI/M^b/c. 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 post-departure 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\times 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^$a,b$/c that we considered includes most models ([1,2,4,5,6,8,9,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^$a,b$/c, and establish a transition probability matrix (t.p.m.) for the system in Section 3. In Section 4, Section 5 and Section 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^$a,b$/c, there are c independent servers, each serving at the rate $\mu$. The customers arrive at the rate $\lambda$ according to a renewal process with an arbitrary inter-arrival time distribution $A\left(t\right)$. 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.${}^{\prime }$s). We consider the system to be in a steady state with the traffic intensity $\rho =\lambda /\left(bc\mu \right)<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^$a,b$/c, 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 $(S_n,n)$, where $n\ge 0$ is the queue-length and $S_n$ is a supplementary flag defined as

$$S_n=\left\{\begin{array}{cc}I\left(k\right),\hfill & \mathrm{if}\mathrm{k}\mathrm{servers}\mathrm{are}\mathrm{idle},1\le \mathrm{k}\le \mathrm{c},0\le \mathrm{n}\le \mathrm{a}-1,\hfill \\ B,\hfill & \mathrm{if}\mathrm{all}\mathrm{servers}\mathrm{are}\mathrm{busy},\mathrm{n}\ge 0.\hfill \end{array}\right.$$

We define the system as busy if all the servers are busy ($S_n=B$), and idle if at least one server is idle ($S_n=I\left(k\right),k$ is the number of idle servers). The queue-length n can be written as $n=qb+n_0, 0\le n_0\le b-1$, where q is the nearest lower non-negative integer of the fraction n/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.

• $\left[l\right|m;t]$ and $\left[l\right|m]$, where $0\le l\le m\le c$, and there are less than a customers waiting in the queue at the beginning of the period, thus $q=0$. Here,

  $$\left[l\right|m;t]=\left(\genfrac{}{}{0pt}{}{m}{l}\right){(1-e^{-\mu t})}^l{\left(e^{-\mu t}\right)}^{m-l}$$

  is the conditional probability that l of m servers complete services during an inter-arrival period of duration t, given that m servers are busy ($c-m$ servers are idle) at the beginning of the period. Moreover, $\left[l\right|m]$ is defined as

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

  (1)
• $\left\{l\right|c;q\}$ is the conditional probability that l of c servers become idle during an inter-arrival period, given that all c servers are busy at the beginning of the period, and q$(q\ge 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\mu$ until time V has elapsed. When all c servers are busy, the number of departed batches follows a Poisson process with a rate $c\mu$. The time V is Erlang-distributed, so it is the sum of q exponential random variables with a rate $c\mu$, implying that the probability density function (p.d.f.) of V is given by

  $$p\left(v\right)=\frac{\left(c\mu \right){\left(c\mu v\right)}^{q-1}{e}^{-c\mu v}}{(q-1)!},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 $\left[l\right|c;t-V]$. Therefore

  $$\left\{l\right|c;q\}={\int }_0^{\infty }{\int }_0^t\left(\genfrac{}{}{0pt}{}{c}{l}\right){(1-e^{-(t-v)\mu })}^l{\left(e^{-(t-v)\mu }\right)}^{c-l}\frac{\left(c\mu \right){\left(c\mu v\right)}^{q-1}{e}^{-c\mu v}}{(q-1)!}dvdA\left(t\right).$$

  (2)
• $\left(l\right|c)$ is the conditional probability that l batches complete service during an inter-arrival 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\mu$:

  $$\left(l\right|c)={\int }_0^{\infty }\frac{e^{-c\mu t}{\left(c\mu t\right)}^l}{l!}dA\left(t\right),l\ge 0.$$

  (3)

Remark 1.

• $\left[0\right|c]=\left(0\right|c)={\int }_0^{\infty }{e}^{-c\mu t}dA\left(t\right)\equiv K_0.$
  Though $\left[0\right|c]$ and $\left(0\right|c)$ give identical results, they have totally different meanings. $\left[0\right|c]$ is for the case when $(c-1)$ servers are busy and $(a-1)$ customers are in queue. After one customer arrives, all the servers become busy without any departures during the inter-arrival time. In this situation, the number of customers in the queue must be zero. Moreover, $\left(0\right|c)$ is for the case that all the servers are already busy before an arrival, and no departures happen during an inter-arrival time. In this situation, the queue-length can be any non-negative number.
• It is easy to prove that $\left(l\right|c)=\{0|c;l\}$.

Let $J_r$ be the system state on the arrival of the rth customer who sees n customers in the queue. The entry of the one-step t.p.m. T from state $(S_i,i)$ to state $(S_j,j)$ is

$$\left[T_{(S_i,i),(S_j,j)}\right]=P(J_{r+1}=(S_j,j)|J_r=(S_i,i)),i\ge 0,j\ge 0,$$

implying that the $(r+1)$th arriving customer sees j customers waiting in the queue with the server state $S_j$, given that the previous rth arriving customer saw i customers waiting in the queue with the server state $S_i$.

The Markov chain (see

Table 1. Submatrix ${\mathbf{T}}_{Idle\to Idle}$.

$(S_n,n)$

$\left(I\right(c),0)$

$\left(I\right(c),1)$

⋯

$\left(I\right(c),a-1)$

$\left(I\right(c-1),0)$

$\left(I\right(c-1),1)$

⋯

$\left(I\right(c-1)$, $a-1)$

⋯

$\left(I\right(2),0)$

$\left(I\right(2),1)$

⋯

$\left(I\right(2),a-1)$

$\left(I\right(1),0)$

$\left(I\right(1),1)\cdots$

$\left(I\right(1),a-1)$

$\left(I\right(c),0)$

$\left[0\right|0]$

$⋮$

⋱

$\left(I\right(c),a-2)$

[0|0]

$\left(I\right(c),a-1)$

[1|1]

[0|1]

$\left(I\right(c-1),0)$

[1|1]

[0|1]

$⋮$

$⋮$

⋱

⋱

$\left(I\right(c-1),a-2)$

[1|1]

[0|1]

$\left(I\right(c-1),a-1)$

[2|2]

[1|2]

$⋮$

$⋮$

⋱

⋱

⋱

$\left(I\right(2),0)$

[(c-2)|(c-2)]

[(c-1)|(c-2)]

[0|(c-2)]

$⋮$

$⋮$

⋱

⋱

⋱

⋱

$\left(I\right(2),a-2)$

[(c-2)|(c-2)]

[[c-3|c-2]

[0|(c-2)]

$\left(I\right(2),a-1)$

[(c-1)|(c-1)]

[(c-2)|c-1)]

[1|(c-1)]

[0|(c-1)]

$\left(I\right(1),0)$

[(c-1)|(c-1)]

[(c-2)|(c-1)]

[1|(c-1)]

[0|(c-1)]

$⋮$

$⋮$

⋱

⋱

⋱

⋱

$\left(I\right(1),a-2)$

[(c-1)|(c-1)]

[(c-2)|(c-1)]

[1|(c-1)]

[0|(c-1)]

$\left(I\right(1),a-1)$

[c|c]

[(c-1)|c]

[2|c]

[1|c]

,

Table 2. Submatrix ${\mathbf{T}}_{Busy\to Idle}$.

$(S_n,n)$

$\left(I\right(c),0)$

$\left(I\right(c),1)$

⋯

$\left(I\right(c),a-1)$

$\left(I\right(c-1),0)$

$\left(I\right(c-1),1)$

⋯

$\left(I\right(c-1)$, $a-1)$

⋯

$\left(I\right(2),0)$

$\left(I\right(2),1)$

⋯

$\left(I\right(2),a-1)$

$\left(I\right(1),0)$

$\left(I\right(1),1)\cdots$

$\left(I\right(1),a-1)$

$(B,0)$

[c|c]

[(c-1)|c]

[2|c]

[1|c]

$⋮$

$⋮$

⋱

⋱

⋱

⋱

$(B,a-2)$

[c|c]

[(c-1)|c]

[2|c]

[1|c]

$(B,a-1)$

{c|c;1}

{(c-1)|c;1}

{2|c;1}

{1|c;1}

$⋮$

$⋮$

⋱

⋱

⋱

⋱

$(B,b-1)$

{c|c;1}

{(c-1)|c;1}

{2|c;1}

{1|c;1}

$(B,b)$

{c|c;1}

{(c-1)|c;1}

{2|c;1}

{1|c;1}

$⋮$

$⋮$

⋱

⋱

⋱

⋱

$(B,b+a-2)$

{c|c;1}

{(c-1)|c;1}

{2|c;1}

{1|c;1}

$(B,b+a-1)$

{c|c;2}

{(c-1)|c;2}

{2|c;2}

{1|c;2}

$⋮$

$⋮$

⋱

⋱

⋱

⋱

$(B,2b-1)$

{c|c;2}

{(c-1)|c;2}

{2|c;2}

{1|c;2}

$(B,2b)$

{c|c;2}

{(c-1)|c;2}

{2|c;2}

{1|c;2}

$⋮$

$⋮$

⋱

⋱

⋱

⋱

$(B,(q-1\left)b\right)$

{c|c;q-1}

{(c-1)|c;q-1}

{2|c;q-1}

{1|c;q-1}

$⋮$

$⋮$

⋱

⋱

⋱

⋱

$(B,(q-1)b+a-2)$

{c|c;q-1}

{(c-1)|c;q-1}

{2|c;q-1}]

{1|c;q-1}

$(B,(q-1)b+a-1)$

{c|c;q}

{(c-1)|c;q}

{2|c;q}

{1|c;q}

$⋮$

$⋮$

⋱

⋱

⋱

⋱

$(B,qb-1)$

{c|c;q}

{(c-1)|c;q}

{2|c;q}

{1|c;q}

$(B,qb)$

{c|c;q}

{(c-1)|c;q}

{2|c;q}

{1|c;q}

$⋮$

$⋮$

⋱

⋱

⋱

⋱

,

Table 3. Submatrix ${\mathbf{T}}_{Busy\to Busy}$.

$(S_n,n)$	$(B,0)$	$(B,1)$	$(B,2)$	⋯	$(B,a)$	⋯	$(B,b)$	⋯
$(B,0)$		(0|c)
$(B,1)$			(0|c)
$⋮$	$⋮$			⋱		⋱		⋱
$(B,a-2)$
$(B,a-1)$	(1|c)				(0|c)
$⋮$	$⋮$			⋱		⋱		⋱
$(B,b-1)$	(1|c)						(0|c)
$(B,b)$		(1|c)
$⋮$	$⋮$			⋱		⋱		⋱
$(B,b+a-2)$
$(B,b+a-1)$	(2|c)
$⋮$	$⋮$			⋱		⋱		⋱
$(B,2b-1)$	(2|c)
$(B,2b)$		(2|c)
$⋮$	$⋮$			⋱		⋱		⋱
$(B,(q-1\left)b\right)$		((q-1)|c)
$⋮$	$⋮$			⋱		⋱		⋱
$(B,(q-1)b+a-2)$
$(B,(q-1)b+a-1)$	(q|c)				((q-1)|c)
$⋮$	$⋮$			⋱		⋱		⋱
$(B,qb-1)$	(q|c)						((q-1)|c)
$(B,qb)$		(q|c)
$⋮$	$⋮$			⋱		⋱		⋱

and

Table 4. Submatrix ${\mathbf{T}}_{Idle\to Busy}$.

$(S_n,n)$	$(B,0)$	$(B,1)$	⋯	$(B,a-1)$	⋯	$(B,b)$	⋯
$\left(I\right(c),0)$
$⋮$
$\left(I\right(c),a-2)$
$\left(I\right(c),a-1)$
$\left(I\right(c-1),0)$
$⋮$
$\left(I\right(c-1),a-2)$
$\left(I\right(c-1),a-1)$
$⋮$
$\left(I\right(2),0)$
$⋮$
$\left(I\right(2),a-2)$
$\left(I\right(2),a-1)$
$\left(I\right(1),0)$
$⋮$
$\left(I\right(1),a-2)$
$\left(I\right(1),a-1)$	[0|c]

) for this system is two-dimensional rather than the usual one-dimensional. The t.p.m. can be formed as four sub-matrices, which are shown in Table 1, Table 2, Table 3 and Table 4.

We describe the four sub-matrices that form the t.p.m.

$$\mathbf{T}=\left[\begin{array}{cc}{\mathbf{T}}_{Idle\to Idle}& {\mathbf{T}}_{Idle\to Busy}\\ {\mathbf{T}}_{Busy\to Idle}& {\mathbf{T}}_{Busy\to Busy}\end{array}\right].$$

(4)

(I)

${\mathbf{T}}_{Idle\to Idle}$. In this situation, the number of customers waiting in queue is less than a. Assume that there are $k_i$ servers idle at the beginning of the inter-arrival time period, and $k_j$ servers idle at the end of the inter-arrival time period, $1\le k_i\le k_j\le c$.

$$[T_{(S_i,i),(S_j,j)}]=\left\{\begin{array}{cc}[T_{(I\left(k_i\right),i),(I\left(k_j\right),i+1)}]=\left[(k_j-k_i)\right|(c-k_i)]\hfill & \mathrm{if}0\le i<a-1,j=i+1,\hfill \\ [T_{(I\left(k_i\right),a-1),(I\left(k_i\right),0)}]=\left[(k_j-k_i+1)\right|(c-k_i+1)]\hfill & \mathrm{if}i=a-1,j=0.\hfill \end{array}\right.$$

(5)

(II)

${\mathbf{T}}_{Busy\to Idle}$. All the servers are busy at the beginning of the period, and $k(1\le k\le c)$ servers are idle at the end of the period, implying that the number of customers in the queue, say j, at the end of the period, must be less than a, i.e., $j<a$. In a manner similar to what we define for $n=qb+n_0,0\le n_0\le b-1$, we need to arrange i customers who are waiting in queue, with FCFS discipline, into q full-size batches and a batch holding the remainders, i.e., $i=qb+i_0,0\le i_0\le b-1$.

$$[T_{(S_i,i),(S_j,j)}]=\left\{\begin{array}{cc}[T_{(B,i),\left(I\right(k),i+1)}]=\left[k\right|c]\hfill & \mathrm{if}0\le i<a-1,j=i+1,\hfill \\ [T_{(B,qb+i_0),(I\left(k\right),i_0+1)}]=\left\{k\right|c;q\}\hfill & \mathrm{if}0\le i_0<a-1,q\ge 1,j=i_0+1,\hfill \\ [T_{(B,qb+i_0),(I\left(k\right),0)}]=\left\{k\right|c;q+1\}\hfill & \mathrm{if}a-1\le i_0\le b-1,q\ge 0,j=0.\hfill \end{array}\right.$$

(6)

(III)

${\mathbf{T}}_{Idle\to Busy}$. The system is idle at the beginning of the time period. After one customer arrives, all the servers become busy and are still busy at the end of the time period. This case appears only if the number of customers waiting in queue is $a-1$, and there is only one server idle at the beginning of the time period.

$$[T_{(S_i,i),(S_j,j)}]=\left\{\begin{array}{cc}[T_{\left(I\right(1),a-1),(B,0)}]=\left[0\right|c]\hfill & \mathrm{if}i=j-1,j=0,\hfill \\ [T_{(I\left(k_i\right),i),(B,j)}]=0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(7)

(IV)

${\mathbf{T}}_{Busy\to Busy}$. All the servers are busy from the beginning to the end of the period, and the number of batches served in time t follows the Poisson process with rate $c\mu$.

$$[T_{(S_i,i),(S_j,j)}]=\left\{\begin{array}{cc}[T_{(B,qb+i_0),(B,(q-l)b+i_0+1)}]=\left(l\right|c)\hfill & \mathrm{if}0\le i_0<b-1,0\le l\le q,\hfill \\ \hfill & i=qb+i_0,j=(q-l)b+1,\hfill \\ [T_{(B,qb+i_0),(B,0)}]=(q+1|c)\hfill & \mathrm{if}a-1\le i_0\le b-1andj=0,q\ge 0.\hfill \end{array}\right.$$

(8)

Finally, $\left[T_{(S_i,i),(S_j,j)}\right]=0$ if $j>i+1$ is true for all of the above I–IV cases. By using identities 1 and 2, it can be easily proven that the sum of all the entries in t.p.m. equals one.

Identity 1.

${\sum }_{l=1}^c\left\{l\right|c;q\}+{\sum }_{i=0}^q\left(i\right|c)=1$ for $q>0.$ This equation shows that the sum of all the conditional probabilities in each row of t.p.m. (when the initial system state is busy) equals one.

Proof. 

$$\begin{array}{cc}\hfill \sum _{l=1}^c\left\{l|c;q\right\}& ={\int }_0^{\infty }{\int }_0^t\sum _{l=1}^c\left(\genfrac{}{}{0pt}{}{c}{l}\right){(1-e^{-(t-v)\mu })}^l{\left(e^{-(t-v)\mu }\right)}^{c-l}\times \frac{\left(c\mu \right){\left(c\mu v\right)}^{q-1}{e}^{-c\mu v}}{(q-1)!}dvdA\left(t\right)\hfill \\ & ={\int }_0^{\infty }{\int }_0^t(1-e^{-c\mu (t-v)})\frac{\left(c\mu \right){\left(c\mu v\right)}^{q-1}{e}^{-c\mu v}}{(q-1)!}dvdA\left(t\right)\hfill \\ & =\underset{Term1}{\underbrace{{\int }_0^{\infty }{\int }_0^t\frac{\left(c\mu \right){\left(c\mu v\right)}^{q-1}{e}^{-c\mu v}}{(q-1)!}dvdA\left(t\right)}}-\underset{Term2}{\underbrace{{\int }_0^{\infty }{\int }_0^t\frac{\left(c\mu \right){\left(c\mu v\right)}^{q-1}{e}^{-c\mu t}}{(q-1)!}dvdA\left(t\right).}}\hfill \end{array}$$

“Term 1” in the above equation can be simplified as $1-{\sum }_{i=0}^{q-1}\left(i\right|c)$ by using the results that the CDF of Erlang is $1-{\sum }_{i=0}^{q-1}\frac{{\left(c\mu t\right)}^i{e}^{-c\mu t}}{i!}$ and $\left(i|c\right)={\int }_0^{\infty }\frac{{\left(c\mu t\right)}^i{e}^{-c\mu t}}{i!}dA\left(t\right).$ “Term 2” can be simplified to $\left(q|c\right)$. Combining these two terms gives ${\sum }_{l=1}^c\left\{l|c;q\right\}=1-\sum _{i=0}^q\left(i\right|c).$ □

Identity 2.

${\sum }_{i=m}^c\left[(i-m)\right|(c-m)]=1,0\le m\le c$. This equation shows that, when the initial system state is idle, the sum of all the conditional probabilities in each row of t.p.m. equals one.

Proof. 

$$\begin{gathered} {\sum }_{i=m}^c\left[(i-m)\right|(c-m)] \\ ={\int }_0^{\infty }{\sum }_{i=m}^c\left(\genfrac{}{}{0pt}{}{c-m}{i-m}\right){(1-e^{-\mu t})}^{i-m}{\left(e^{-\mu t}\right)}^{c-i}dA\left(t\right) \\ ={\int }_0^{\infty }dA\left(t\right)=1. \end{gathered}$$

□

Since the Markov chain under consideration is irreducible, positive recurrent and aperiodic, it has a limiting distribution if and only if $\rho =\lambda /bc\mu <1$. In view of this, $\underset{r\to \infty }{lim}P(J_r=(S_n,n))=X(S_n,n)$ exists. In this case, the limiting distribution is given by $\mathbf{X}=\mathbf{XT}$ where $\mathbf{T}$ is t.p.m. defined in (4), and the vector $\mathbf{X}$ has the form

$$\begin{array}{cc}\hfill \mathbf{X}=\left[X\right(I\left(c\right),0),\cdots ,X(I\left(c\right),a-1),\cdots ,& X\left(I\right(1),0),\cdots ,\hfill \\ \hfill & X\left(I\right(1),a-1),X(B,0),\cdots X(B,1),\cdots ],\hfill \end{array}$$

(9)

where $X\left(I\right(k),n),0\le n<a$ and $X(B,n),n\ge 0$, respectively, denote the p.a.e. unnormalized probabilities that an arriving customer sees n customers in queue, k of c servers idle, and n customers in queue, with all servers busy. If such a vector $\mathbf{X}$ exists, it will be the vector of the steady state p.a.e. probabilities up to some normalizing constant.

4. Queue-Length Distributions at Pre-Arrival Epoch

4.1. The Busy Server Probabilities

When all the servers are busy during an inter-arrival time period, for the queueing model GI/M^$a,b$/c, the service times for batches are i.i.d.r.v.${}^{\prime }$ s, having exponential distributions. Thus, the number of batches that complete service during an arbitrary inter-arrival time will have a Poisson distribution, which implies that the probability of l service completions during an inter-arrival time A is $\left(l\right|c)$, and the probability generating function (p.g.f.) of $\left(l\right|c)$ is

$$D\left(z\right)=\sum _{l=0}^{\infty }\left(l\right|c)z^l=\overline{a}\left(c\mu (1-z)\right),$$

(10)

where $\overline{a}\left(\alpha \right)$ is the Laplace–Stieltjes transform (L.-S.T.) of $A\left(t\right)$, i.e., $\overline{a}\left(\alpha \right)={\int }_0^{\infty }exp(-\alpha t)dA\left(t\right)$ and

$$K_0=\overline{a}\left(c\mu \right)={\int }_0^{\infty }exp(-c\mu t)dA\left(t\right).$$

(11)

Theorem 1.

For the queueing system GI/M^$a,b$/c, in the steady state case, the busy-server probabilities of queue length at pre-arrival epoch are given by $P^{-}(B,n)=X(B,n)/C_N=w^n/C_N,n\ge 0$, where w is a real root inside the unit circle of equation $D\left(z^b\right)=z=\overline{a}\left(c\mu (1-z^b)\right)$ and $C_N$ is a normalizing constant given by $C_N={\sum }_{j=1}^c{\sum }_{i=0}^{a-1}X\left(I\right(j),i)+\frac{1}{1-w.}$

Proof. 

When the system is busy and n customers are waiting in the queue, it is evident from t.p.m. that

$$X(B,n)=\sum _{j=0}^{\infty }\left(j\right|c\left)X\right(B,jb+n-1),n>0.$$

(12)

To solve the difference Equation (12), in the same manner as by Chaudhry and Madill [5], a solution of the form $X(B,n)=z^n(z\ne 0),n\ge 1$ is assumed. For more details on this method, one may see Chaudhry and Templeton ( [14], page 350). Substituting $X(B,n)=z^n$ into Equation (12), we have

$$z^n=\sum _{j=0}^{\infty }{z}^{jb+n-1}\left(j\right|c)=z^{n-1}\sum _{j=0}^{\infty }\left(j\right|c)z^{jb}=z^{n-1}D\left(z^b\right).$$

(13)

Combining this with Equation (10), and simplifying, we obtain the root equation

$$D\left(z^b\right)=z=\overline{a}\left(c\mu (1-z^b)\right).$$

(14)

By Rouché’s theorem, it can be shown that Equation (14) has a real root w inside the unit circle if $\rho =\frac{\lambda }{bc\mu }<1$. Once the root w is found, $X(B,1),X(B,2),\cdots$ can be obtained by using $X(B,n)=w^n,n\ge 1$.

Next, we can solve for $X(B,0)$ from Equation (12) by setting $n=1$,

$$X(B,1)=w=\sum _{j=0}^{\infty }\left(j\right|c\left)X\right(B,jb),$$

implying

$$w=\left(0\right|c)X(B,0)+\left(1\right|c)w^b+\left(2\right|c)w^{2b}+\cdots =D\left(w^b\right).$$

(15)

Combining (10), (14), and (15), we conclude $X(B,0)=1$. This implies that the assumption $X(B,n)=z^n$ is true even for $n=0.$

Finally, $P^{-}(B,n)$ can be obtained as the normalized $X(B,n)$ by dividing a normalizing constant $C_N$ (see Equations (23) and (24)). □

4.2. The Idle Server Probabilities

The idle server unnormalized probabilities $X\left(I\right(c),0),\cdots ,X\left(I\right(c),a-1),\cdots ,X\left(I\right(1),0),\cdots ,X\left(I\right(1),a-1)$ can be obtained by $c\times a$ linear equations generated from the t.p.m. In fact, there are $‘‘c\times a+1$” equations, with (as usual) one being redundant.

These “$c\times a+1$” equations are

$$\begin{array}{cc}\hfill X(B,0)=& X(I\left(1\right),a-1)\left[0\right|c]+\sum _{i=1}^{\infty }\left(i\right|c)\sum _{l=a}^{b}X(B,(i-1)b+l-1),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill X\left(I\right(k),j)=& \sum _{m=1}^{k}X\left(I\right(m),j-1)\left[(k-m)\left|\right(c-m)\right]+X(B,j-1)\left[k\right|c]+\hfill \\ \hfill & \sum _{i=1}^{\infty }X(B,ib+j-1)\left\{k\right|c;i\},\hfill \end{array}$$

(17)

$$\begin{array}{c}\hfill \mathrm{and}\\ \hfill X\left(I\right(k),0)=& \sum _{m=1}^{k+1}X\left(I\right(m),a-1)\left[\right(k-m+1\left)\right|(c-m+1)]+\hfill \\ \hfill & \sum _{i=1}^{\infty }\left\{k\right|c;i\}\sum _{l=a}^{b}X(B,(i-1)b+l-1),\hfill \end{array}$$

(18)

where $1\le j\le a-1,1\le k\le c$ and $X\left(I\right(c+1),a-1)=0$.

Remark 2.

The $c\times a$ idle server unknown probabilities (unnormalized)

$$\begin{array}{c}\hfill \left[X\right(I\left(c\right),0),\cdots ,X(I\left(c\right),a-1),\cdots ,X(I\left(1\right),0),\cdots ,X(I\left(1\right),a-1\left)\right]\end{array}$$

can be obtained simultaneously by using the above $c\times a$ equations. However, large values of c or a may cause a computational problem, since the last terms in both Equations (17) and (18) are infinite series related to complex double integrals $\left\{k\right|c;i\}$ (defined in Equation (2)). In general, when we operate on an infinite series without a closed form, the series has to be truncated. Therefore, the result is approximated as we lose the tails due to this truncation. To fix these problems, we want to simplify Equations (17) and (18) by deriving a closed form for these series. Before we move on, we need to prove the following two lemmas.

Lemma 1.

The last term in Equation (16)

$$\begin{array}{cc}\hfill & \sum _{i=1}^{\infty }\left(i\right|c)\sum _{l=a}^{b}X(B,(i-1)b+l-1)=\frac{w^{a-b-1}-1}{1-w}(w-K_0).\hfill \end{array}$$

Proof. 

$$\begin{array}{cc}\hfill & \sum _{i=1}^{\infty }\left(i\right|c)\sum _{l=a}^{b}X(B,(i-1)b+l-1)\hfill \\ \hfill & =w^{a-b-1}\sum _{i=1}^{\infty }\left(i\right|c)\sum _{k=0}^{b-a}{w}^{ib+k}=w^{a-b-1}\frac{1-w^{b-a+1}}{1-w}\sum _{i=1}^{\infty }\left(i\right|c)w^{ib}\hfill \\ \hfill & =w^{a-b-1}\frac{1-w^{b-a+1}}{1-w}(D\left(w^b\right)-K_0)=\frac{w^{a-b-1}-1}{1-w}(w-K_0)\hfill \end{array}$$

by using $\left(0\right|c)=K_0$, and Equation (13). □

Lemma 2.

Define $J\left(k\right)={\sum }_{i=1}^{\infty }{w}^{ib}\left\{k\right|c;i\},$ and

$$J\left(k\right)=c\mu w^b{\int }_0^{\infty }{\int }_0^t\left(\genfrac{}{}{0pt}{}{c}{k}\right){(1-e^{-(t-v)\mu })}^k{\left(e^{-(t-v)\mu }\right)}^{c-k}{e}^{-c\mu v(1-w^b)}dvdA\left(t\right).$$

(19)

Proof. 

$$\begin{array}{cc}& \sum _{i=1}^{\infty }{w}^{ib}\left\{k\right|c;i\}\hfill \\ & =\sum _{i=1}^{\infty }{w}^{ib}{\int }_0^{\infty }{\int }_0^t\left(\genfrac{}{}{0pt}{}{c}{k}\right){(1-e^{-(t-v)\mu })}^k{\left(e^{-(t-v)\mu }\right)}^{c-k}\frac{\left(c\mu \right){\left(c\mu v\right)}^{i-1}{e}^{-c\mu v}}{(i-1)!}dvdA\left(t\right)\hfill \\ & ={\int }_0^{\infty }{\int }_0^t\left(\genfrac{}{}{0pt}{}{c}{k}\right){(1-e^{-(t-v)\mu })}^k{\left(e^{-(t-v)\mu }\right)}^{c-k}\sum _{i=1}^{\infty }{w}^{ib}\frac{\left(c\mu \right){\left(c\mu v\right)}^{i-1}{e}^{-c\mu v}}{(i-1)!}dvdA\left(t\right)\hfill \\ & =c\mu w^b{\int }_0^{\infty }{\int }_0^t\left(\genfrac{}{}{0pt}{}{c}{k}\right){(1-e^{-(t-v)\mu })}^k{\left(e^{-(t-v)\mu }\right)}^{c-k}{e}^{-c\mu v(1-w^b)}\underset{=1, \mathrm{Poisson}\mathrm{p}.\mathrm{m}.\mathrm{f}}{\underbrace{\sum _{i=1}^{\infty }\frac{{\left(c\mu v{w}^b\right)}^{i-1}{e}^{-c\mu v{w}^b}}{(i-1)!}}}dvdA\left(t\right).\hfill \end{array}$$

□

Theorem 2.

For the queueing system GI/M^$a,b$/c, in the steady state case, the idle server probabilities of queue length at the pre-arrival epoch are given by $P^{-}(I\left(k\right),n)=X(I\left(k\right),n)/C_N,0\le n<a-1,1\le k\le c$, where $C_N$ is a normalizing constant given by $C_N={\sum }_{j=1}^c{\sum }_{i=0}^{a-1}X\left(I\right(j),i)+\frac{1}{1-w}$ and $X\left(I\right(k),n)$ satisfy the following equations:

$$\left(i\right)X(I\left(1\right),a-1)=\frac{1}{(1-w)K_0}(1-w^{a-b}+K_0{w}^{a-b-1}-K_0),$$

(20)

$$\left(ii\right)X(I\left(k\right),j)=\sum _{m=1}^{k}X\left(I\right(m),j-1)\left[\right(k-m\left)\right|(c-m)]+w^{j-1}(\left[k\right|c]+J\left(k\right)),1<j<a-1,$$

(21)

$$\left(iii\right)X(I\left(k\right),0)=\sum _{m=1}^{k+1}X\left(I\right(m),a-1)\left[\right(k-m+1\left)\right|(c-m+1)]+\frac{w^{a-b-1}-1}{1-w}J\left(k\right).$$

(22)

Proof. 

(i) Using Lemma 1 and $\left[0\right|c]=K_0$, we can rewrite Equation (16) and directly solve for $X\left(I\right(1),a-1)$.

(ii) and (iii) Using Theorem 1, replacing $X(B,j-1)$ by $w^{j-1}$, $X(B,ib+j-1)$ by $w^{ib+j-1}$, and $X(B,(i-1)b+l-1)$ by $w^{(i-1)b+l-1}$, then applying the result of Lemma 2, we can rewrite Equations (17) and (18) as Equations (21) and (22), respectively.

We first solved $X\left(I\right(1),a-1)$ using Equation (20), and then solved other idle server probabilities recursively by using Equations (21) and (22). For more details on this, see the algorithm developed in Appendix A. □

Finally, we obtained all queue-length probabilities, and needed to normalize the vector

$$\begin{array}{cc}\hfill \mathbf{X}=\left[X\right(I\left(c\right),0),\cdots ,X(I\left(c\right),a-1),\cdots ,& X\left(I\right(1),0),\cdots ,\hfill \\ & X\left(I\right(1),a-1),X(B,0),\cdots X(B,1),\cdots ].\hfill \end{array}$$

by dividing a normalizing constant $C_N$, which is given by

$$C_N=\sum _{j=1}^c\sum _{i=0}^{a-1}X\left(I\right(j),i)+\sum _{i=0}^{\infty }X(B,i)=\sum _{j=1}^c\sum _{i=0}^{a-1}X\left(I\right(j),i)+\frac{1}{1-w}.$$

(23)

Define ${\mathbf{P}}^{-}$ as the vector of normalized p.a.e. such that

$${\mathbf{P}}^{-}=\frac{\mathbf{X}}{C_N}.$$

(24)

Further, $P^{-}(I\left(k\right),n)$ and $P^{-}(B,n)$, respectively, are normalized p.a.e. probabilities and represent that k of the c servers are idle, $0\le n<a-1$, and all servers are busy, $n\ge 0$.

4.3. Special Cases

4.3.1. Single-Server Probabilities for GI/M^$a,b$/1

The system GI/M^$a,b$/1 is a special case of GI/M^$a,b$/c when $c=1$.

(A)

When $c=1$, the root Equation (14) is simplified to $D\left(z\right)=z=\overline{a}\left(\mu (1-z)\right)$, which agrees with the root equation in the work by Chaudhry and Madill [5]; consequently, the same results of $X(B,0),\cdots ,X(B,1),\cdots ,X(B,M)$ can be obtained.

(B)

Moreover, $k=m=c=1,\left[0\right|0]=1,\left[1\right|1]=1-\left[0\right|1]=1-K_0$, and ${\sum }_{i=1}^{\infty }{w}^{ib}\left\{1\right|1;i\}=\frac{1}{(1-w^b)}\left(w^b-w+(1-w^b)K_0\right)$. Equation (21) can be simplified to

$$\begin{array}{cc}\hfill X\left(I\right(1),j)& =X(I\left(1\right),j-1)+w^{j-1}(1-K_0+\frac{1}{(1-w^b)}(w^b-w+(1-w^b)K_0))\hfill \\ & =X(I\left(1\right),j-1)+w^{j-1}\frac{1-w}{1-w^b}.\hfill \end{array}$$

This agrees with the equation in Chaudhry and Madill [5] for solving the idle server probabilities.

4.3.2. Multi-Server Queueing System GI/M^b/c

The system GI/M^b/c is a special case of GI/M^$a,b$/c when $a=1$.

In GI/M^b/c, the system is idle only if there is no customer waiting in queue. Instead of evaluating the queue-length distributions, Chaudhry and Templeton [14] consider the distribution for the number of customers in the system for GI/M^b/c without considering the server being busy or idle. The numerical results for the system GI/M^b/c are also not available. We can see that our model includes this model as a special case, it not only produces the numerical solutions for the queue-length distributions, but also the information of the server utilization.

5. Queue-Length Distributions at Random Epoch

We are now interested in knowing the probability that the system will be in a given state at a random epoch (r.e.) in time. A random epoch is said to occur at the end of a random period of time, R, since the last p.a.e. From renewal theory, the probability associated with R, $dR\left(t\right)$ is given by $dR\left(t\right)=\lambda (1-A(t\left)\right)dt,t>0$ (see Chaudhry and Templeton [14]). Proceeding in a manner directly analogous to that used for developing $\left(l\right|c)$, $\left[l\right|m]$ and $\left\{l\right|c;q\}$, where the services are considered during the inter-arrival time A (see Equations (1)–(3)), ${\left(l\right|c)}_R$, ${\left[l\right|m]}_R$ and ${\left\{l\right|c;q\}}_R$ are defined as the probabilities that such services take place during time R. The p.g.f. of ${\left(l\right|c)}_R$ (see proof in Appendix B) is

$$D_R\left(z\right)=\sum _{l=0}^{\infty }{\left(l\right|c)}_R{z}^l=\frac{\rho b}{1-z}[1-\overline{a}\left(c\mu (1-z)\right)],$$

(25)

and

$${\left(0\right|c)}_R={\left[0\right|c]}_R=\lambda {\int }_0^{\infty }exp(-c\mu t)(1-A(t\left)\right)dt=\rho b(1-K_0).$$

(26)

Similar to the definition for the p.a.e probability vector ${\mathbf{P}}^{-}$ in Equation (24), we define $\mathbf{P}$ as the vector of the r.e. probabilities, such that

$$\mathbf{P}=\left[P\right(I\left(c\right),0),\cdots ,P(I\left(c\right),a-1),\cdots ,P(I\left(1\right),0),\cdots ,P(I\left(1\right),a-1),P(B,0),\cdots P(B,1),\cdots ],$$

where $P\left(I\right(k),n),0\le n<a$ and $P(B,n),n\ge 0$, respectively, denote the r.e. probabilities that, at the end of a random period of time R after arrival, k of the c servers are idle, $0\le n<a-1$ customers are in the queue, and all servers are busy, $n\ge 0$ customers are in the queue. The forms of the t.p.m. $\mathbf{T}$ in Table 1, Table 2, Table 3 and Table 4 contain all of the information required on transitions within the queueing system in a period measured from the last p.a.e. The nature of the entries in the t.p.m. serve to indicate the probabilities associated with the transitions. Thus, if the limiting distribution is ${\mathbf{P}}^{-}={\mathbf{P}}^{-}\mathbf{T}$ when the timeframe is the inter-arrival time, A, instead of the entries $\left(l\right|c)$, $\left[l\right|m]$ and $\left\{l\right|c;q\}$, the entries ${\left(l\right|c)}_R$, ${\left[l\right|m]}_R$ and ${\left\{l\right|c;q\}}_R$ are used with the timeframe, R, and $\mathbf{P}={\mathbf{P}}^{-}{\mathbf{T}}_{\mathrm{R}}$, where the newly formed t.p.m. ${\mathbf{T}}_{\mathrm{R}}$ describes how the steady-state p.a.e. probabilities are transformed into steady-state probabilities for the system at a random epoch after the last p.a.e.

Remark 3.

Similar to those in the p.a.e. systems, it can be proven that the following three equations still hold for the case of r.e. systems:

• ${\left[0\right|c]}_R={\left(0\right|c)}_R\equiv \rho b(1-K_0)$ (see Equation (26));
• ${\sum }_{l=1}^c{\left\{l\right|c;q\}}_R+{\sum }_{i=0}^q{\left(i\right|c)}_R=1$ for $q>0$; and
• ${\sum }_{i=m}^c{\left[(i-m)\right|(c-m)]}_R=1,0\le m\le c$.

Thus, the sum of entries in each row of t.p.m. ${\mathbf{T}}_{\mathrm{R}}$ equals one.

5.1. The Busy-Server Probabilities

The busy-server r.e. probabilities $P(B,n),n\ge 0$ can be calculated in a similar manner as the queue-length distributions at the pre-arrival epoch described in Section 4.1. Here, we derive the closed-form busy-server probability distribution of the queue length at a random epoch. The probabilities $P(B,n),n\ge 0$ can be obtained using Equations (27) and (28) (see below). Since both are in terms of the root w, the calculations become extremely simple. The key idea to derive these two equations is based on the relations between two probabilities: $P(B,n)$ and $P^{-}(B,n),n\ge 0$.

Theorem 3.

For the queueing system GI/M^$a,b$/c, in the steady state case, the busy-server probabilities of the queue length at the random epoch are given by

(i) 

$P(B,n)=\frac{1}{C_N}\frac{\rho b(1-w)w^{n-1}}{1-w^b},n>0.$

(ii) 

$P(B,0)=\frac{\rho b(1-K_0)}{C_N(1-w)K_0}(1-w^{a-b})+\frac{\rho b(w^{a-b-1}-1)}{C_N(1-w^b)}.$

Proof. 

(i) At the end of a random period of time R after arrival, if all servers are busy and the waiting line is not empty ($n>0$), then the sizes for those batches that were taken into service during time R must be maximum ($=b$, full batch size). Since the queue length at a pre-arrival epoch will be $n-1+mb,m\ge 0$, it leads to r.e. probabilities as

$$P(B,n)=\sum _{m=0}^{\infty }{\left(m\right|c)}_R{P}^{-}(B,mb+n-1),n>0.$$

By using the fact that $P^{-}(B,n)=w^n,$ and Equations (14) and (25), we have

$$\begin{array}{cc}\hfill P(B,n)& =\frac{1}{C_N}\sum _{m=0}^{\infty }{\left(m\right|c)}_R{w}^{mb+n-1}\hfill \\ \hfill & =\frac{1}{C_N}\frac{\rho b(1-w)w^{n-1}}{1-w^b},n>0.\hfill \end{array}$$

(27)

(ii) In this situation, the queue length is empty at a random time while all the servers are busy, then the size for the last batch into service can be any number between $[a,b]$, and the servers at the moment when the last customer arrives are either all busy or one idle. Combining all of these possibilities, using Equation (20) and the following equation

$$\begin{array}{cc}\hfill \sum _{i=1}^{\infty }{\left(i\right|c)}_R\sum _{j=a}^b{P}^{-}\left(B,(i-1)b+j-1\right)& =\frac{w^{a-b-1}-1}{C_N(1-w)}\left[\sum _{i=0}^{\infty }{\left(i\right|c)}_R{w}^{ib}-{\left(0\right|c)}_R\right]\hfill \\ & =\frac{w^{a-b-1}-1}{C_N(1-w)}\left(\frac{\rho b(1-w)}{1-w^b}-\rho b(1-K_0)\right),\hfill \end{array}$$

$P(B,0)$ can be expressed as

$$\begin{array}{cc}\hfill & P(B,0)\hfill \\ \hfill & ={\left(0\right|c)}_R{P}^{-}(I\left(1\right),a-1)+\sum _{i=1}^{\infty }{\left(i\right|c)}_R\sum _{j=a}^b{P}^{-}(B,(i-1)b+j-1)\hfill \\ \hfill & =\frac{\rho b(1-K_0)}{C_N(1-w)K_0}(1-w^{a-b}+K_0{w}^{a-b-1}-K_0)+\frac{w^{a-b-1}-1}{C_N(1-w)}(\frac{\rho b(1-w)}{1-w^b}-\rho b(1-K_0))\hfill \\ \hfill & =\frac{\rho b(1-K_0)}{C_N(1-w)K_0}(1-w^{a-b})+\frac{\rho b(w^{a-b-1}-1)}{C_N(1-w^b)}.\hfill \end{array}$$

(28)

□

At the end of a random period of time R after arrival, if all servers are busy, the queue length n $(n\ge 0)$ distribution can be evaluated by using Equations (27) and (28). In this case, both the results are in closed-form in terms of the root w.

5.2. The Idle Server Probabilities

Corollary 1.

The idle server r.e. probabilities $P\left(I\right(k),n),0\le n<a$ can be obtained by using Theorem 2. The Equations (30) and (31) (see below) are modified from Equations (21) and (22) in Theorem 2 by replacing the term $\left[l|m\right]$ with ${\left[l|m\right]}_R$, and normalizing the probabilities from $X\left(I\right(m),j-1)$ to $P^{-}(I\left(m\right),j-1),1<j<a$. Moreover, we redefine $J_R\left(k\right)$ as

$$\begin{array}{cc}\hfill J_R\left(k\right)& =\sum _{i=1}^{\infty }{w}^{ib}{\left\{k|c;i\right\}}_R\hfill \\ \hfill & =c\lambda \mu w^b\left({\int }_0^{\infty }{\int }_0^t\left(\genfrac{}{}{0pt}{}{c}{k}\right){(1-e^{-(t-v)\mu })}^k{\left(e^{-(t-v)\mu }\right)}^{c-k}{e}^{-c\mu v(1-w^b)}(1-A\left(t\right))dvdt\right)\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill Then,& P(I\left(k\right),0)=\sum _{m=1}^{k+1}{P}^{-}(I\left(m\right),a-1){\left[k-m+1|c-m+1\right]}_R+\frac{w^{a-b-1}-1}{1-w}{J}_R\left(k\right),\hfill \end{array}$$

(30)

$$P(I\left(k\right),j)=\sum _{m=1}^k{P}^{-}(I\left(m\right),j-1){\left[k-m|c-m\right]}_R+w^{j-1}\left({\left[k|c\right]}_R+J_R\left(k\right)\right),$$

(31)

where $1\le j\le a-1,1\le k\le c,P^{-}(I(c+1),a-1)=0$.

5.3. The Special Case: E_$\eta$/M^$a,b$/c Queue

The system E_$\eta$/M^$a,b$/1 is a special case of GI/M^$a,b$/c when the inter-arrival time is Erlang (with $\eta$ phase)-distributed. Then the root Equation (14) can be simplified to

$${\left(\frac{\eta \rho b}{\eta \rho b+1-z^b}\right)}^{\eta }-z=0.$$

By replacing $dA\left(t\right)$ with $\frac{{\left(\lambda \eta \right)}^{\eta }{t}^{\eta -1}{e}^{-\lambda \eta t}}{(\eta -1)!}dt$, we can calculate p.a.e. probability distributions for both busy and idle servers by using the algorithm introduced in Appendix A. Then the r.e. probability distributions can be obtained by using Equations (27)–(31).

Sim [10] solved the $\eta$-phase Erlangian arrivals system E_$\eta$/M^$a,b$/c only for the probabilities at r.e. and discussed the results in the context of transportation systems. Our algorithms can not only solve the systems with general inter-arrival time distributions, but also provide the solutions at different epochs. Our numerical results agree with those provided by Sim [10].

6. Queue-Length Distributions at Post-Departure Epoch

In this section, we derive the probabilities for the state of the system immediately after a real service completion takes place. It was assumed that no time elapsed after the server completed a batch before accepting a quorum-complete batch from the queue. Thus, the post-departure epoch (p.d.e.) occurred immediately after a server had either reduced the queue or became idle.

To find the p.d.e. probabilities, we need to first define an epoch—a pre-service completion epoch (p.s.e.), i.e., the instant in the time immediately before a real departure occurs (before a real service completes). Then, $P^{S_{-}}(I\left(k\right),n)$ and $P^{S_{-}}(B,n),n\ge 0,1\le k\le c$, respectively, are defined as the probabilities at p.s.e., when there are n customers in queue, k of c servers idle, and n customers in queue, all servers busy. It is apparent that $P^{S_{-}}(I\left(c\right),n)=0$ for any n.

Similarly, we define $P^{+}(I\left(k\right),n),0\le n<a,1\le k\le c$ and $P^{+}(B,n),n\ge 0$, as the probabilities of the queue length at a p.d.e.

Conjecture 1.

The following relationships between p.d.e. and p.s.e. probabilities apply

$$\begin{array}{cc}\hfill & P^{+}(I\left(k\right),n)=P^{S_{-}}(I(k-1),n),0\le n\le a-1,2\le k\le c\hfill \\ \hfill & P^{+}(I\left(1\right),n)=P^{S_{-}}(B,n),0\le n\le a-1,\hfill \end{array}$$

(32)

and

$$\begin{array}{cc}\hfill & P^{+}(B,n)=P^{S_{-}}(B,n+b),n\ge 1,\hfill \\ \hfill & P^{+}(B,0)=\sum _{n=a}^b{P}^{S_{-}}(B,n).\hfill \end{array}$$

(33)

Corollary 2.

$P^{S_{-}}(I\left(k\right),n),0\le n<a,1\le k\le c$ and $P^{S_{-}}(B,n),n\ge 0$ satisfy the following equations:

$$\begin{array}{cc}\hfill & P^{S_{-}}(I\left(k\right),n)=\frac{P\left(I\right(k),n)}{1-{\sum }_{i=0}^{a-1}P\left(I\right(c),i)},0\le n\le a-1,1\le k\le c-1,\hfill \\ \hfill & P^{S_{-}}(B,n)=\frac{P(B,n)}{1-{\sum }_{i=0}^{a-1}P\left(I\right(c),i)},n\ge 0.\hfill \end{array}$$

(34)

Proof. 

When the service time distribution is exponential, service completions, real or potential, occur at random epochs. The probabilities, $P^{S_{-}}(I\left(k\right),n),0\le n<a,1\le k\le c$ and $P^{S_{-}}(B,n),n\ge 0$ can be found by conditioning the r.e. probabilities to ensure that at least one server is busy. Thus, using the results of r.e. probabilities given in Theorem 3, we can obtain p.d.e. probabilities for both busy and idle servers from Equations (32)–(34). □

Remark 4.

(i) 

When we set $c=1$, these probabilities agree with those of Chaudhry and Madill [5] for the system GI/M^$a,b$/1.

(ii) 

As a check on the algebra, also useful as a computational check, we note that, using (32)–(34),

$$\begin{array}{cc}\hfill \sum _{k=1}^c\sum _{j=0}^{a-1}{P}^{+}(I\left(k\right),j)+\sum _{j=0}^{\infty }{P}^{+}(B,j)& =\frac{{\sum }_{k=1}^{c-1}{\sum }_{j=0}^{a-1}P(I\left(k\right),j)+{\sum }_{j=0}^{\infty }P(B,j)}{1-{\sum }_{i=0}^{a-1}P\left(I\right(c),i)}\hfill \\ & =\frac{1-{\sum }_{i=0}^{a-1}P\left(I\right(c),i)}{1-{\sum }_{i=0}^{a-1}P\left(I\right(c),i)}=1,\hfill \end{array}$$

as it should be.

7. Numerical Results

In this section, we present some numerical results for various inter-arrival time distributions such as $\eta$-phase Erlang (E_$\eta$), deterministic (D), and uniform (U). All the examples we considered have the same mean value of the inter-arrival time $E\left(A\right)=1/\lambda$. The root equation (see Equation (14)), probability density functions (p.d.f.) of inter-arrival time A, and p.d.f. of a random period time R for these three distributions are summarized in

Table 5. Root Equations, p.d.f.s of $A\left(t\right)$, $R\left(t\right)$, and mean value of of $A\left(t\right)$ for E_$\eta$/M^$a,b$/c, D/M^$a,b$/c and U/M^$a,b$/c.

Inter-arrival time distributions	Root Equations (Equation (14))	p.d.f. of $A\left(t\right)$	p.d.f. of $R\left(t\right)$	$E\left(A\right)$
$\eta$-phase Erlang	${\left(\frac{\eta \rho b}{\eta \rho b+1-z^b}\right)}^{\eta }-z=0$	$\frac{{\left(\lambda \eta \right)}^{\eta }{t}^{\eta -1}exp(-\lambda \eta t)}{(\eta -1)!}$	$\lambda {\sum }_{n=0}^{\eta -1}\frac{{\left(\lambda \eta t\right)}^{n}exp(-\lambda \eta t)}{n!}$	$1/\lambda$
Deterministic	$exp(-\frac{1-z^b}{\rho b})-z=0$	$\delta (t-1/\lambda )$	$\left\{\begin{array}{c}\lambda ,ift<\frac{1}{\lambda }\hfill \\ 0,ift\ge \frac{1}{\lambda }\hfill \end{array}\right.$	$1/\lambda$
Uniform	$\frac{exp(-\frac{1-z^b}{\rho b})}{\phi c\mu (1-z^b)}\times [exp(\phi c\mu (1-z^b)/2)-exp(-\phi c\mu (1-z^b))/2]-z=0,$$\phi =t_2-t_1,\mathrm{is}\mathrm{the}\mathrm{interval}\mathrm{width}$	$1/\phi$	$\left\{\begin{array}{c}\lambda ,ift<t_1\hfill \\ \frac{1}{\phi }+\frac{\lambda }{2}-\frac{\lambda t}{\phi },if{t}_1\le t<t_2\hfill \\ 0,ift\ge t_2\hfill \end{array}\right.$$t_1=\frac{1}{\lambda }-\frac{\phi }{2},t_2=\frac{1}{\lambda }+\frac{\phi }{2}$	$1/\lambda$

.

Besides the calculations for the queue-length probabilities at the pre-arrival, random, and post-departure epochs for both idle and busy systems, we also considered the performance measures, such as the mean (denoted as LQ^e) and the standard deviations (denoted as SDLQ^e) of the queue length; the mean (denoted as $E^e\left[I\left(k\right)\right]$) and variance (denoted as $Va{r}^e\left[I\left(k\right)\right]$ ) of the idle servers. The symbol “e” denotes the epoch state, which can be pre-arrival (e = “−”), random (e = “ ”), or post-departure (e = “+”). We define $P{B}^e={\sum }_{n=0}^{\infty }{P}^e(B,n)$ as the probability that an arriving customer sees the system busy at e epoch, and $P{I}^e={\sum }_{n=0}^{a-1}{\sum }_{k=1}^c{P}^e(I\left(k\right),n)$ is the probability that the system is idle at e epoch. The probabilities of the queue length at three different epochs are presented in closed form. Since most of these probabilities are irrational, for computational purposes, we need to set the precision $\epsilon$. Throughout all computations in the following examples, we use $\epsilon$ = 10${}^{-20}$ as the precision. Due to the rounding error, the sum of the probabilities may not be one.

The results of the E₆/M^5,10/5 queue with traffic intensities $\rho =0.1,0.5,0.9$ for both busy and idle servers at pre-arrival epoch are presented in

Table 6. Distribution of queue lengths at pre-arrival epochs for the busy system E₆/M^5,10/5, $\rho =0.1,0.5,0.9,\epsilon ={10}^{-20}$.

n	$P^{-}(B,n)$
	$\rho =0.1$	$\rho =0.5$	$\rho =0.9$
0	1.12 × 10${}^{-5}$	0.0715625	0.0198544
1	4.45 × 10${}^{-6}$	0.0612334	0.0194400
2	1.76 × 10${}^{-6}$	0.0523951	0.0190343
3	7.00 × 10${}^{-7}$	0.0448325	0.0186371
4	2.77 × 10${}^{-7}$	0.0383615	0.0182481
5	1.10 × 10${}^{-7}$	0.0328244	0.0178673
⋮	⋮	⋮	⋮
10	1.08 × 10${}^{-9}$	0.015056	0.0160790
⋮	⋮	⋮	⋮
50	9.28 × 10${}^{-26}$	2.95 × 10${}^{-5}$	0.0069163
⋮	⋮	⋮	⋮
296	⋮	6.55 × 10${}^{-22}$	3.86 × 10${}^{-5}$
⋮	⋮	⋮	⋮
2184	⋮	⋮	1.96 × 10${}^{-22}$
$P{B}^{-}$	0.0000186	0.4957989	0.9513294
$P{I}^{-}$	0.9999815	0.5042011	0.0486706
SUM	1.0000001	1.0000000	1.0000000

and

Table 7. Distribution of queue lengths at the pre-arrival epochs for the idle system E₆/M^5,10/5, $\rho =0.1,0.5,0.9,\epsilon ={10}^{-20}$.

$\rho =0.9$		n					Probabilities of k servers idle
		0	1	2	3	4
$I\left(k\right)$	1	0.0045133	0.0062023	0.0077054	0.0090392	0.0102190	0.0376792
	2	0.0008471	0.0012652	0.0017914	0.0024039	0.0030832	0.0093908
	3	0.0000996	0.0001656	0.0002602	0.0003883	0.0005533	0.0014670
	4	0.0000058	0.0000116	0.0000210	0.0000351	0.0000554	0.0001289
	5	0.0000001	30.0000003	0.0000007	0.0000013	0.0000023	0.0000047
# in queue		0.0054659	0.0076450	0.0097787	0.0118678	0.0139132	SUM: 0.0486706
$\rho =0.5$		n					Probabilities of k servers idle
		0	1	2	3	4
$I\left(k\right)$	1	0.0476528	0.0526108	0.0551242	0.0557986	0.0551149	0.2663013
	2	0.0231853	0.0281785	0.0331282	0.0377295	0.0417906	0.1640121
	3	0.0064926	0.0089895	0.0118683	0.0150732	0.0185253	0.0609489
	4	0.0008326	0.0014196	0.0021965	0.0031834	0.0043942	0.0120263
	5	0.0000258	0.0000722	0.0001465	0.0002567	0.0004115	0.0009126
# in queue		0.0781891	0.0912706	0.1024636	0.1120414	0.1202365	SUM: 0.5042011
$\rho =0.1$		n					Probabilities of k servers idle
		0	1	2	3	4
$I\left(k\right)$	1	0.0005343	0.0002563	0.0001226	0.0000585	0.0000279	0.0009995
	2	0.009807	0.0057356	0.0033333	0.0019273	0.0011097	0.0219129
	3	0.0620315	0.0455535	0.0329312	0.0235189	0.0166373	0.1806724
	4	0.1067216	0.1061743	0.1006805	0.0923661	0.0827349	0.4886774
	5	0.0208943	0.0422757	0.0629306	0.0821285	0.0994901	0.0009996
# in queue		0.1999887	0.1999954	0.1999982	0.1999993	0.1999999	SUM: 0.9999815

, respectively. When we set the number of servers to 1, our results match with those obtained for E₆/M^5,10/1 by Chaudhry et al. [5].

We considered three systems E₆/M^a,10/5, D/M^a,10/5, and U/M^a,10/5 ($t_1=0.875/\lambda ,t_2=1.125/\lambda ,\phi =0.25/\lambda )$. All three systems have the same mean value of inter-arrival time $E\left(A\right)=1/\lambda$. In

Table 8. Comparison of performance measures of E₆/M^a,10/5, D/M^a,10/5, and U/M^a,10/5 for idle servers, a = 1, 4, 7, $\rho$ = 0.1, 0.5, 0.9, $\epsilon$ = 10${}^{-20}$.

	$\rho$	a	$P{I}^{-}$	$L_Q^{-}$	$SDL{Q}^{-}$	$E^{-}\left[I\left(k\right)\right]$	$PI$	$L_Q$	$SDLQ$	$E\left[I\right(k\left)\right]$	$P{I}^{+}$	$L_Q^{+}$	$SDL{Q}^{+}$	$E^{+}\left[I\left(k\right)\right]$
E6/Ma, 10/5	0.1	1	0.5918	0.2683	0.7409	1.0643	0.4824	0.4082	0.8821	0.8093	0.7533	0.0000	0.0095	1.5576
		4	0.9998	1.5000	1.1181	3.8502	0.9997	1.5000	1.1181	3.7500	1.0000	1.3564	1.1081	4.4644
		7	1.0000	3.0000	2.0000	4.3430	1.0000	3.0000	2.0000	4.2857	1.0000	2.3522	1.9248	4.8083
	0.5	1	0.0205	5.8067	6.3981	0.0223	0.0107	6.2025	6.4242	0.0115	0.1047	1.3049	3.8825	0.1163
		4	0.3396	4.4994	5.6122	0.4972	0.3140	4.7317	5.7188	0.4523	0.6218	1.8572	3.1865	1.0733
		7	0.7796	3.7588	3.6688	1.6043	0.7670	3.8177	3.7469	1.5536	0.9205	3.0814	2.5287	2.4557
	0.9	1	0.0015	46.8450	47.4120	0.0016	0.0007	47.2590	47.4150	0.0008	0.0137	38.2720	46.5250	0.0145
		4	0.0285	45.6330	47.3290	0.0340	0.0249	46.0310	47.3480	0.0295	0.0991	37.4040	46.1810	0.1287
		7	0.1149	41.9250	46.7130	0.1636	0.1093	42.2760	46.7680	0.1543	0.2295	34.6680	44.9320	0.3837
D/Ma, 10/5	0.1	1	0.6002	0.2327	0.6704	1.0461	0.4601	0.3998	0.8399	0.7354	0.7471	0.0000	0.0063	1.4800
		4	0.9999	1.5000	1.1180	3.8708	0.9998	1.5000	1.1800	3.7500	1.0000	1.3556	1.1077	4.4753
		7	1.0000	3.0000	2.0000	4.3548	1.0000	3.0000	2.0000	4.2857	1.0000	2.3415	1.9220	4.8142
	0.5	1	0.0182	5.6832	6.2591	0.0194	0.0076	6.1610	6.2870	0.0080	0.0924	1.2518	3.7652	0.1004
		4	0.3390	4.4119	5.4834	0.4914	0.3076	4.6909	5.6106	0.4374	0.6211	1.8256	3.0937	1.0579
		7	0.7842	3.7155	3.5716	1.6102	0.7690	3.7841	3.6631	1.5491	0.9229	3.0714	2.4851	2.4549
	0.9	1	0.0013	45.0490	46.6060	0.0014	0.0005	46.5470	46.6100	0.0005	0.0120	37.5580	45.7220	0.0125
		4	0.0280	44.8720	46.5260	0.0331	0.0237	45.3500	46.5470	0.0278	0.0983	36.7210	45.3830	0.1261
		7	0.1151	41.2030	45.9100	0.1630	0.1083	41.6240	45.9770	0.1518	0.2302	34.0240	44.1370	0.3819
U/Ma, 10/5	0.1	1	0.6000	0.2338	0.6726	1.0468	0.4608	0.4001	0.8412	0.7378	0.7473	0.0000	0.0064	1.4825
		4	0.9999	1.5000	1.1180	3.8702	0.9998	1.5000	1.1180	3.7500	1.0000	1.3557	1.1077	4.4799
		7	1.0000	3.0000	2.0000	4.3544	1.0000	3.0000	2.0000	4.2857	1.0000	2.3419	1.9221	4.8140
	0.5	1	0.0182	5.6870	6.2634	0.0195	0.0076	6.1623	6.2913	0.0081	0.0928	1.2535	3.7689	0.1001
		4	0.3390	4.4117	5.4874	0.4916	0.3079	4.6921	5.6140	0.4378	0.6211	1.8266	3.0966	1.0584
		7	0.7840	3.7168	3.5746	1.6100	0.7690	3.7852	3.6657	1.5492	0.9229	3.0717	2.4865	2.4549
	0.9	1	0.0013	46.0740	46.6320	0.0014	0.0005	46.5690	46.6350	0.0005	0.0120	37.5810	45.7480	0.0125
		4	0.0280	44.8960	46.5500	0.0332	0.0238	45.3720	46.5720	0.0279	0.0983	36.7420	45.4080	0.1262
		7	0.1151	41.2260	45.9350	0.1631	0.1083	41.6440	46.0020	0.1519	0.2302	34.0450	44.1620	0.3819

, we present the performance measures for these three systems for idle servers at three different epochs with varied $a=1,4,7$ and $\rho =0.1,0.5,0.9$. In Figure 1, we compare the performance of D/M^4,10/5 for busy servers at pre-arrival epochs with $\rho =0.1,0.3,0.5,0.7$ and $0.9$. In Figure 2, we compare the performance of U/M^a,10/5 for busy servers at pre-arrival epochs with $a=1,4,7$.

8. Conclusions

The queue GI/M^$a,b$/c was successfully investigated by using the two-dimensional embedded Markov chain. Simple and exact analyses to determine queue-length distributions are presented. An algorithm was derived for the analysis of the steady state behaviour of the system. Our recursive solution approach is not only very efficient, but also accurate by providing the exact queue-length probabilities at p.a.e. In a similar manner, we studied the queue-length distribution at r.e. and derived closed-form formulae in terms of the root w for evaluating the exact queue-length probabilities at r.e. We also obtained the probabilities of p.d.e. through the relations between r.e. and p.d.e. The results for this system were provided numerically by considering three inter-arrival time distributions—Erlang, deterministic, and uniform. The work on higher order moments and other distributions can be conducted similarly.

There are two special features in this work. The first is the effort to express the important results in closed form; the second is the development of the methodology and algorithms to efficiently derive accurate results. The models under consideration were validated by using MAPLE to obtain numerical results with sufficient accuracy and trivial computational costs.