**Remark 1.**


Let *Jr* be the system state on the arrival of the *r*th customer who sees *n* customers in the queue. The entry of the one-step t.p.m. *T* from state (*Si*, *i*) to state (*Sj*, *j*) is

$$[T\_{(S\_i,i),(S\_{\backslash i},j)}] = P(f\_{r+1} = (S\_{j\prime}j)|f\_r = (S\_{i\prime}i)), \qquad i \ge 0, j \ge 0, i$$

implying that the (*r* + 1)th arriving customer sees *j* customers waiting in the queue with the server state *Sj*, given that the previous *r*th arriving customer saw *i* customers waiting in the queue with the server state *Si*.

The Markov chain (see Tables 1–4) 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 Tables 1–4.

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

$$\mathbf{T} = \begin{bmatrix} \mathbf{T}\_{Idllc \to Idlc} & \mathbf{T}\_{Idllc \to Busy} \\ \mathbf{T}\_{Busy \to Idlc} & \mathbf{T}\_{Busy \to Busy} \end{bmatrix} \tag{4}$$

(I) **T***Idle*→*Idle*. In this situation, the number of customers waiting in queue is less than *a*. Assume that there are *ki* servers idle at the beginning of the inter-arrival time period, and *kj* servers idle at the end of the inter-arrival time period, 1 ≤ *ki* ≤ *kj* ≤ *c*.

$$[T\_{(S\_{\bar{i}},i),(S\_{\bar{j}},\bar{i})}] = \begin{cases} [T\_{(l(k\_i),i),(l(k\_j),i+1)}] = [(k\_j - k\_i)|(c - k\_i)] & \text{if } 0 \le i < a - 1, j = i + 1, \\ [T\_{(l(k\_i),a - 1),(l(k\_i),0)}] = [(k\_j - k\_i + 1)|(c - k\_i + 1)] & \text{if } i = a - 1, j = 0. \end{cases} \tag{5}$$

(II) **<sup>T</sup>***Busy*→*Idle*. All the servers are busy at the beginning of the period, and *k*(1 ≤ *k* ≤ *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 ≤ *n*0 ≤ *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 ≤ *i*0 ≤ *b* − 1.

$$\begin{cases} \begin{aligned} \left[T\_{(S\_{\boldsymbol{\nu}}\boldsymbol{l}),(S\_{\boldsymbol{\rho}}\boldsymbol{l})}\right] = \begin{cases} \left[T\_{(B,\boldsymbol{l}),(I(k),\boldsymbol{l}+1)}\right] = \left[k|\boldsymbol{c}\right] & \text{if } 0 \le \boldsymbol{i} < a-1, \boldsymbol{j} = \boldsymbol{i}+1, \\\left[T\_{(B,qh+\boldsymbol{l}\_{0}),(I(k),\boldsymbol{j}+1)}\right] = \left\{k|\boldsymbol{c};q\right\} & \text{if } 0 \le \boldsymbol{i}\_{0} < a-1, q \ge 1, \boldsymbol{j} = \boldsymbol{i}\_{0}+1, \\\left[T\_{(B,qh+\boldsymbol{l}\_{0}),(I(k),\boldsymbol{\mathcal{O}})}\right] = \left\{k|\boldsymbol{c};q+1\right\} & \text{if } a-1 \le \boldsymbol{i}\_{0} \le b-1, q \ge 0, \boldsymbol{j}=0. \end{aligned} \end{cases} \end{cases}$$

(III) **<sup>T</sup>***Idle*→*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\_{(\mathcal{S}\_i,i),(\mathcal{S}\_j,j)}] = \begin{cases} [T\_{(I(1),\mu-1),(\mathcal{B},0)}] = [0|\mathcal{c}] & \text{if } i=j-1, j=0, \\ [T\_{(I(k\_i),i),(\mathcal{B},j)}] = 0 & \text{otherwise.} \end{cases} \tag{7}$$

(IV) **<sup>T</sup>***Busy*→*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 *<sup>c</sup>μ*.

$$[T\_{(S\_l j),(S\_{\tilde{j}})}] = \begin{cases} [T\_{(B\_\eta \mathfrak{b} + i\_0),(B\_\eta (q-l)b + i\_0 + 1)}] = (l|c) & \text{if } 0 \le i\_0 < b - 1, 0 \le l \le q\_\nu \\ & i = qb + i\_0, j = (q-l)b + 1, \\\ [T\_{(B\_\eta \mathfrak{b} + i\_0),(B\_\eta 0)}] = (q+1|c) & \text{if } a - 1 \le i\_0 \le b - 1 \text{ and } j = 0, q \ge 0. \end{cases} \tag{8}$$

Finally, [*<sup>T</sup>*(*Si*,*<sup>i</sup>*),(*Sj*,*<sup>j</sup>*)] = 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.** ∑*cl*=<sup>1</sup>{*l*|*<sup>c</sup>*; *q*} + <sup>∑</sup>*qi*=<sup>0</sup>(*i*|*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{split} \sum\_{l=1}^{\varepsilon} \{l|c;q\} &= \int\_{0}^{\infty} \int\_{0}^{t} \sum\_{l=1}^{\varepsilon} \binom{\varepsilon}{l} (1-e^{-(t-\upsilon)\mu})^{l} (e^{-(t-\upsilon)\mu})^{\varepsilon-l} \times \frac{(c\mu)(c\mu\upsilon)^{q-1}e^{-c\mu\upsilon}}{(q-1)!} d\upsilon dA(t) \\ &= \int\_{0}^{\infty} \int\_{0}^{t} (1-e^{-\varepsilon\mu(t-\upsilon)}) \frac{(c\mu)(c\mu\upsilon)^{q-1}e^{-c\mu\upsilon}}{(q-1)!} d\upsilon dA(t) \\ &= \underbrace{\int\_{0}^{\infty} \int\_{0}^{t} \frac{(c\mu)(c\mu\upsilon)^{q-1}e^{-c\mu\upsilon}}{(q-1)!} d\upsilon dA(t)}\_{\text{Term 1}} - \underbrace{\int\_{0}^{\infty} \int\_{0}^{t} \frac{(c\mu)(c\mu\upsilon)^{q-1}e^{-c\mu t}}{(q-1)!} d\upsilon dA(t)}\_{\text{Term 2}}. \end{split}$$

"Term 1" in the above equation can be simplified as 1 − ∑*<sup>q</sup>*−<sup>1</sup> *i*=0 (*i*|*c*) by using the results that the CDF of Erlang is 1 − ∑*<sup>q</sup>*−<sup>1</sup> *i*=0 (*cμ<sup>t</sup>*)*ie*−*<sup>c</sup>μ<sup>t</sup> i*! and (*i*|*c*) = ∞0 (*cμ<sup>t</sup>*)*ie*−*<sup>c</sup>μ<sup>t</sup> i* ! *dA*(*t*). "Term 2" can be simplified to (*q*|*c*). Combining these two terms gives *c* ∑ *l*=1 {*l*|*c*; *q*} = 1 − *q* ∑ *i*=0 (*i*|*c*).

**Identity 2.** ∑*ci*=*<sup>m</sup>*[(*<sup>i</sup>* − *m*)|(*c* − *m*)] = 1, 0 ≤ *m* ≤ *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.*

$$\begin{array}{l} \textbf{Proof.} \ \sum\_{i=m}^{\zeta} |(i-m)|(c-m) | \\ = \int\_{0}^{\infty} \sum\_{i=m}^{\zeta} \binom{c-m}{i-m} (1-e^{-\mu t})^{i-m} (e^{-\mu t})^{c-i} dA(t) \\ = \int\_{0}^{\infty} dA(t) = 1. \end{array}$$






**Table 3.** Submatrix **<sup>T</sup>***Busy*→*Busy*.

Since the Markov chain under consideration is irreducible, positive recurrent and aperiodic, it has a limiting distribution if and only if *ρ* = *λ*/*bcμ* < 1. In view of this, lim *r*→∞ *P*(*Jr* = (*Sn*, *n*)) = *<sup>X</sup>*(*Sn*, *n*) exists. In this case, the limiting distribution is given by **X** = **XT** where **T** is t.p.m. defined in (4), and the vector **X** has the form

$$\begin{aligned} \mathbf{X} = \left[ X(I(\mathbf{c}), 0), \dots, X(I(\mathbf{c}), a-1), \dots, X(I(1), 0), \dots \right] \\ X(I(1), a-1), X(B, 0), \dots, X(B, 1), \dots \right], \end{aligned} \tag{9}$$

where *<sup>X</sup>*(*I*(*k*), *<sup>n</sup>*), 0 ≤ *n* < *a* and *<sup>X</sup>*(*<sup>B</sup>*, *<sup>n</sup>*), *n* ≥ 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 **X** exists, it will be the vector of the steady state p.a.e. probabilities up to some normalizing constant.


**Table 4.** Submatrix **<sup>T</sup>***Idle*→*Busy*.

#### **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*<sup>a</sup>*, *<sup>b</sup>*/*<sup>c</sup>*, the service times for batches are i.i.d.r.v.- 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 (*l*|*c*), and the probability generating function (p.g.f.) of (*l*|*c*) is

$$D(z) = \sum\_{l=0}^{\infty} (l|c)z^l = \mathfrak{a}(c\mu(1-z)),\tag{10}$$

where *a*¯(*α*) is the Laplace–Stieltjes transform (L.-S.T.) of *<sup>A</sup>*(*t*), i.e., *a*¯(*α*) = ∞0 exp(−*α<sup>t</sup>*)*dA*(*t*) and 

$$\mathbb{K}\_0 = \mathbb{d}(c\mu) = \int\_0^\infty \exp(-c\mu t) dA(t). \tag{11}$$

**Theorem 1.** *For the queueing system GI/Ma*, *b/c, in the steady state case, the busy-server probabilities of queue length at pre-arrival epoch are given by <sup>P</sup>*−(*<sup>B</sup>*, *n*) = *<sup>X</sup>*(*<sup>B</sup>*, *n*)/*CN* = *<sup>w</sup><sup>n</sup>*/*CN*, *n* ≥ 0*, where w is a real root inside the unit circle of equation <sup>D</sup>*(*z<sup>b</sup>*) = *z* = *a*¯(*cμ*(<sup>1</sup> − *z<sup>b</sup>*)) *and CN is a normalizing constant given by CN* = ∑*cj*=<sup>1</sup> ∑*<sup>a</sup>*−<sup>1</sup> *i*=0 *<sup>X</sup>*(*I*(*j*), *i*) + 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|c\rangle X(B, jb+n-1), \qquad n>0. \tag{12}$$

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

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

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

$$D(z^b) = z = \mathbb{1}(c\mu(1 - z^b)).\tag{14}$$

By Rouché's theorem, it can be shown that Equation (14) has a real root *w* inside the unit circle if *ρ* = *λbcμ* < 1. Once the root *w* is found, *<sup>X</sup>*(*<sup>B</sup>*, <sup>1</sup>), *<sup>X</sup>*(*<sup>B</sup>*, <sup>2</sup>), ··· can be obtained by using *<sup>X</sup>*(*<sup>B</sup>*, *n*) = *<sup>w</sup>n*, *n* ≥ 1.

Next, we can solve for *<sup>X</sup>*(*<sup>B</sup>*, 0) from Equation (12) by setting *n* = 1,

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

implying

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

Combining (10), (14), and (15), we conclude *<sup>X</sup>*(*<sup>B</sup>*, 0) = 1. This implies that the assumption *<sup>X</sup>*(*<sup>B</sup>*, *n*) = *z<sup>n</sup>* is true even for *n* = 0.

Finally, *<sup>P</sup>*−(*<sup>B</sup>*, *n*) can be obtained as the normalized *<sup>X</sup>*(*<sup>B</sup>*, *n*) by dividing a normalizing constant *CN* (see Equations (23) and (24)).

#### *4.2. The Idle Server Probabilities*

The idle server unnormalized probabilities *<sup>X</sup>*(*I*(*c*), <sup>0</sup>), ··· , *<sup>X</sup>*(*I*(*c*), *a* − <sup>1</sup>), ··· , *<sup>X</sup>*(*I*(1), <sup>0</sup>), ··· , *<sup>X</sup>*(*I*(1), *a* − 1) can be obtained by *c* × *a* linear equations generated from the t.p.m. In fact, there are "*c* × *a* + 1" equations, with (as usual) one being redundant.

These "*c* × *a* + 1" equations are

$$X(B,0) = X(I(1), a-1)[0|c] + \sum\_{i=1}^{\infty} \left(i|c|\right) \sum\_{l=a}^{b} X(B, (i-1)b+l-1),\tag{16}$$

$$X(I(k),j) = \sum\_{m=1}^{k} X(I(m),j-1)[(k-m)|(c-m)] + X(B,j-1)[k|c] + \sum\_{m=1}^{\infty} X(B,ib+j-1)\{k|c;i\},\tag{17}$$

$$\sum\_{i=1}^{\infty} X(B,ib+j-1)\{k|c;i\},\tag{17}$$

and

$$X(I(k),0) = \sum\_{m=1}^{k+1} X(I(m), a-1)[(k-m+1)(c-m+1)]+$$

$$\sum\_{i=1}^{\infty} \{k \vert c; i\} \sum\_{l=a}^{b} X(B\_r(i-1)b+l-1),\tag{18}$$

where 1 ≤ *j* ≤ *a* − 1, 1 ≤ *k* ≤ *c* and *<sup>X</sup>*(*I*(*c* + <sup>1</sup>), *a* − 1) = 0.

**Remark 2.** *The c* × *a idle server unknown probabilities (unnormalized)*

$$\left[ \left[ X(I(c),0), \dots, X(I(c),a-1), \dots, X(I(1),0), \dots, X(I(1),a-1) \right] \right]$$

*can be obtained simultaneously by using the above c* × *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* {*k*|*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)

$$\sum\_{i=1}^{\infty} \left( i \vert c \right) \sum\_{l=a}^{b} X(B\_{\prime}(i-1)b + l - 1) = \frac{w^{a-b-1} - 1}{1 - w} (w - K\_0).$$

**Proof.**

$$\begin{aligned} &\sum\_{i=1}^{\infty} \left( i|c| \right) \sum\_{l=a}^{b} X(B\_r(i-1)b+l-1) \\ &= w^{a-b-1} \sum\_{i=1}^{\infty} \left( i|c| \right) \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|c| \right) w^{ib} \\ &= w^{a-b-1} \frac{1-w^{b-a+1}}{1-w} (D(w^b) - K\_0) = \frac{w^{a-b-1}-1}{1-w} (w-K\_0) \end{aligned}$$

by using (0|*c*) = *K*0, and Equation (13).

**Lemma 2.** *Define J*(*k*) = ∑∞*<sup>i</sup>*=<sup>1</sup> *wib*{*k*|*c*; *i*}, *and*

$$J(k) = c\mu w^b \int\_0^\infty \int\_0^t \binom{c}{k} (1 - e^{-(t-v)\mu})^k (e^{-(t-v)\mu})^{c-k} e^{-c\mu v (1 - w^b)} dv dA(t). \tag{19}$$

**Proof.**

$$\begin{split} &\sum\_{l=1}^{\infty} w^{lb} \{k \vert c; i\} \\ &= \sum\_{l=1}^{\infty} w^{lb} \int\_{0}^{\infty} \int\_{0}^{t} \binom{c}{k} (1 - e^{-(t-\nu)\mu})^{k} (e^{-(t-\nu)\mu})^{c-k} \frac{(c\mu)(c\mu\upsilon)^{i-1} \varepsilon^{-c\mu\upsilon}}{(i-1)!} d\upsilon dA(t) \\ &= \int\_{0}^{\infty} \int\_{0}^{t} \binom{c}{k} (1 - e^{-(t-\nu)\mu})^{k} (e^{-(t-\nu)\mu})^{c-k} \sum\_{l=1}^{\infty} w^{lb} \frac{(c\mu)(c\mu\upsilon)^{i-1} \varepsilon^{-c\mu\upsilon}}{(i-1)!} d\upsilon dA(t) \\ &= c\mu w^{b} \int\_{0}^{\infty} \int\_{0}^{t} \binom{c}{k} (1 - e^{-(t-\nu)\mu})^{k} (e^{-(t-\nu)\mu})^{c-k} e^{-c\mu\upsilon(1-\mu^{\flat})} \underbrace{\sum\_{l=1}^{\infty} \frac{(c\mu\upsilon\nu^{b})^{l-1} \varepsilon^{-c\mu\upsilon\nu^{b}}}{(i-1)!}}\_{=: \text{\$\!-1\$, \text{\$\!-1\$}\$\!-\!\!$$

**Theorem 2.** *For the queueing system GI/Ma*, *b/c, in the steady state case, the idle server probabilities of queue length at the pre-arrival epoch are given by <sup>P</sup>*−(*I*(*k*), *n*) = *<sup>X</sup>*(*I*(*k*), *<sup>n</sup>*)/*CN*, 0 ≤ *n* < *a* − 1, 1 ≤ *k* ≤ *c, where CN is a normalizing constant given by CN* = ∑*cj*=<sup>1</sup> ∑*<sup>a</sup>*−<sup>1</sup> *i*=0 *<sup>X</sup>*(*I*(*j*), *i*) + 1 1−*<sup>w</sup> and <sup>X</sup>*(*I*(*k*), *n*) *satisfy the following equations:*

$$(i)\ X(I(1),a-1) = \frac{1}{(1-w)K\_0}(1-w^{a-b}+K\_0w^{a-b-1}-K\_0),\tag{20}$$

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

$$\text{(iii)}\quad X(I(k),0) = \sum\_{m=1}^{k+1} X(I(m),a-1)[(k-m+1)(c-m+1)] + \frac{w^{a-b-1}-1}{1-w}I(k). \tag{22}$$

**Proof.** (i) Using Lemma 1 and [0|*c*] = *K*0, we can rewrite Equation (16) and directly solve for *<sup>X</sup>*(*I*(1), *a* − <sup>1</sup>).

(ii) and (iii) Using Theorem 1, replacing *<sup>X</sup>*(*<sup>B</sup>*, *j* − 1) by *wj*−1, *<sup>X</sup>*(*<sup>B</sup>*, *ib* + *j* − 1) by *wib*+*j*−1, and *<sup>X</sup>*(*<sup>B</sup>*,(*<sup>i</sup>* − 1)*b* + *l* − 1) by *w*(*<sup>i</sup>*−<sup>1</sup>)*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 *<sup>X</sup>*(*I*(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{aligned} \mathbf{X} &= \left[ X(I(\varepsilon), 0), \dots, X(I(\varepsilon), a-1), \dots, X(I(1), 0), \dots \right] \\ &\qquad X(I(1), a-1), X(B, 0), \dots, X(B, 1), \dots \right]. \end{aligned}$$

by dividing a normalizing constant *CN*, which is given by

$$\mathbb{C}\_{N} = \sum\_{j=1}^{c} \sum\_{i=0}^{a-1} X(I(j), i) + \sum\_{i=0}^{\infty} X(B, i) = \sum\_{j=1}^{c} \sum\_{i=0}^{a-1} X(I(j), i) + \frac{1}{1 - w}. \tag{23}$$

Define **P**− as the vector of normalized p.a.e. such that

$$\mathbf{P}^{-} = \frac{\mathbf{x}}{\mathbb{C}\_{N}}.\tag{24}$$

Further, *<sup>P</sup>*−(*I*(*k*), *n*) and *<sup>P</sup>*−(*<sup>B</sup>*, *<sup>n</sup>*), respectively, are normalized p.a.e. probabilities and represent that *k* of the *c* servers are idle, 0 ≤ *n* < *a* − 1, and all servers are busy, *n* ≥ 0.

## *4.3. Special Cases*

4.3.1. Single-Server Probabilities for GI/M*<sup>a</sup>*, *<sup>b</sup>*/1

> The system GI/M*<sup>a</sup>*, *<sup>b</sup>*/1 is a special case of GI/M*<sup>a</sup>*, *<sup>b</sup>*/*c* when *c* = 1.


$$\begin{aligned} X(I(1),j) &= X(I(1),j-1) + w^{j-1}(1 - K\_0 + \frac{1}{(1 - w^b)}(w^b - w + (1 - w^b)K\_0)), \\ &= X(I(1),j-1) + w^{j-1} \frac{1 - w}{1 - w^b}. \end{aligned}$$

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*<sup>b</sup>*/*c*

The system GI/M*<sup>b</sup>*/*c* is a special case of GI/M*<sup>a</sup>*, *<sup>b</sup>*/*c* when *a* = 1.

In GI/M*<sup>b</sup>*/*<sup>c</sup>*, 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*<sup>b</sup>*/*c* without considering the server being busy or idle. The numerical results for the system GI/M*<sup>b</sup>*/*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*(*t*) is given by *dR*(*t*) = *λ*(1 − *<sup>A</sup>*(*t*))*dt*, *t* > 0 (see Chaudhry and Templeton [14]). Proceeding in a manner directly analogous to that used for developing (*l*|*c*), [*l*|*m*] and {*l*|*c*; *q*}, where the services are considered during the inter-arrival time *A* (see Equations (1)–(3)), (*l*|*c*)*<sup>R</sup>*, [*l*|*m*]*R* and {*l*|*c*; *q*}*R* are defined as the probabilities that such services take place during time *R*. The p.g.f. of (*l*|*c*)*R* (see proof in Appendix B) is

$$D\_R(z) = \sum\_{l=0}^{\infty} (l|c)\_R z^l = \frac{\rho b}{1-z} [1 - d(c\mu(1-z))],\tag{25}$$

and

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

Similar to the definition for the p.a.e probability vector **P**− in Equation (24), we define **P** as the vector of the r.e. probabilities, such that

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

where *<sup>P</sup>*(*I*(*k*), *<sup>n</sup>*), 0 ≤ *n* < *a* and *<sup>P</sup>*(*<sup>B</sup>*, *<sup>n</sup>*), *n* ≥ 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 ≤ *n* < *a* − 1 customers are in the queue , and all servers are busy, *n* ≥ 0 customers are in the queue. The forms of the t.p.m. **T** in Tables 1–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 **P**− = **P**−**T** when the timeframe is the inter-arrival time, *A*, instead of the entries (*l*|*c*), [*l*|*m*] and {*l*|*c*; *q*}, the entries (*l*|*c*)*<sup>R</sup>*, [*l*|*m*]*R* and {*l*|*c*; *q*}*R* are used with the timeframe, *R*, and **P** = **P**−**TR**, where the newly formed t.p.m. **TR** 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:*


*Thus, the sum of entries in each row of t.p.m.* **TR** *equals one.*

#### *5.1. The Busy-Server Probabilities*

The busy-server r.e. probabilities *<sup>P</sup>*(*<sup>B</sup>*, *<sup>n</sup>*), *n* ≥ 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 *<sup>P</sup>*(*<sup>B</sup>*, *<sup>n</sup>*), *n* ≥ 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: *<sup>P</sup>*(*<sup>B</sup>*, *n*) and *<sup>P</sup>*−(*<sup>B</sup>*, *<sup>n</sup>*), *n* ≥ 0.

**Theorem 3.** *For the queueing system GI/Ma*, *b/c, in the steady state case, the busy-server probabilities of the queue length at the random epoch are given by*

>

0.

$$P(i) \quad P(B\_\prime n) = \frac{1}{C\_N} \frac{\rho b (1 - w) w^{n-1}}{1 - w^k}, \qquad n \ge 1$$

 *(ii) <sup>P</sup>*(*<sup>B</sup>*, 0) = *ρb*(<sup>1</sup>−*K*0) *CN*(<sup>1</sup>−*<sup>w</sup>*)*<sup>K</sup>*<sup>0</sup> (1 − *wa*−*<sup>b</sup>*) + *<sup>ρ</sup>b*(*wa*−*b*−1−<sup>1</sup>) *CN*(<sup>1</sup>−*w<sup>b</sup>*) .

**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* ≥ 0, it leads to r.e. probabilities as

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

By using the fact that *<sup>P</sup>*−(*<sup>B</sup>*, *n*) = *<sup>w</sup>n*, and Equations (14) and (25), we have

$$\begin{split} P(B, n) &= \frac{1}{\mathbb{C}\_N} \sum\_{m=0}^{\infty} \left( m|c| \mathfrak{c} \right)\_R w^{mb+n-1} \\ &= \frac{1}{\mathbb{C}\_N} \frac{\rho b (1-w) w^{n-1}}{1-w^b}, \qquad n > 0. \end{split} \tag{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{aligned} \sum\_{i=1}^{\infty} (i|c)\_R \sum\_{j=a}^b P^- (B\_\prime(i-1)b + j - 1) &= \frac{w^{a-b-1} - 1}{\mathbb{C}\_N(1-w)} \left[ \sum\_{i=0}^\infty (i|c)\_R w^{ib} - (0|c)\_R \right] \\ &= \frac{w^{a-b-1} - 1}{\mathbb{C}\_N(1-w)} \left( \frac{\rho b (1-w)}{1 - w^b} - \rho b (1 - \mathbb{K}\_0) \right). \end{aligned}$$

*<sup>P</sup>*(*<sup>B</sup>*, 0) can be expressed as

$$\begin{split} &P(B,0) \\ &= (0|c)\_R P^-(I(1), a-1) + \sum\_{l=1}^{\infty} (i|c)\_R \sum\_{j=a}^b P^-(B, (i-1)b+j-1) \\ &= \frac{\rho b(1-K\_0)}{\mathbb{C}\_N(1-w)K\_0} (1-w^{a-b} + K\_0 w^{a-b-1} - K\_0) + \frac{w^{a-b-1}-1}{\mathbb{C}\_N(1-w)} (\frac{\rho b(1-w)}{1-w^b} - \rho b(1-K\_0)) \\ &= \frac{\rho b(1-K\_0)}{\mathbb{C}\_N(1-w)K\_0} (1-w^{a-b}) + \frac{\rho b(w^{a-b-1}-1)}{\mathbb{C}\_N(1-w^b)}. \end{split} \tag{28}$$

At the end of a random period of time *R* after arrival, if all servers are busy, the queue length n (*n* ≥ 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 <sup>P</sup>*(*I*(*k*), *<sup>n</sup>*), 0 ≤ *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*[*l*|*m*] *with* [*l*|*m*]*R, and normalizing the probabilities from <sup>X</sup>*(*I*(*m*), *j* − 1) *to <sup>P</sup>*−(*I*(*m*), *j* − <sup>1</sup>), 1 < *j* < *a. Moreover, we redefine JR*(*k*) *as*

$$\begin{split} f\_{\mathbb{R}}(k) &= \sum\_{l=1}^{\infty} w^{lb} \{ k | c; l \}\_{R} \\ &= c\lambda \mu w^{b} \left( \int\_{0}^{\infty} \int\_{0}^{t} \binom{c}{k} (1 - e^{-(t-\overline{v})\mu})^{k} (e^{-(t-\overline{v})\mu})^{v-k} e^{-\varepsilon \mu v (1 - w^{b})} (1 - A(t)) d\upsilon dt \right) \end{split} \tag{29}$$

$$\text{Then, }\ P(I(k),0) = \sum\_{m=1}^{k+1} P^-(I(m),a-1)[k-m+1|c-m+1]\_R + \frac{\nu^{a-b-1}-1}{1-\nu}f\_R(k),\tag{30}$$

$$P(I(k),j) = \sum\_{m=1}^{k} P^-(I(m),j-1)[k-m|c-m]\_R + w^{j-1}([k|c]\_R + f\_R(k)),\tag{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: <sup>E</sup>η/Ma*, *b/c Queue*

The system <sup>E</sup>*η*/M*<sup>a</sup>*, *<sup>b</sup>*/1 is a special case of GI/M*<sup>a</sup>*, *<sup>b</sup>*/*c* when the inter-arrival time is Erlang (with *η* 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*(*t*) with (*λη*)*η t <sup>η</sup>*−1*e*<sup>−</sup>*λη<sup>t</sup>* (*η*−<sup>1</sup>)! *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 *η*-phase Erlangian arrivals system <sup>E</sup>*η*/M*<sup>a</sup>*, *<sup>b</sup>*/*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, *PS*− (*I*(*k*), *n*) and *PS*− (*<sup>B</sup>*, *<sup>n</sup>*), *n* ≥ 0, 1 ≤ *k* ≤ *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 *PS*− (*I*(*c*), *n*) = 0 for any *n*.

Similarly, we define *<sup>P</sup>*+(*I*(*k*), *<sup>n</sup>*), 0 ≤ *n* < *a*, 1 ≤ *k* ≤ *c* and *<sup>P</sup>*+(*<sup>B</sup>*, *<sup>n</sup>*), *n* ≥ 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{aligned} P^+(I(k), n) &= P^{S\_-}(I(k-1), n), \quad 0 \le n \le a-1, 2 \le k \le c\\ P^+(I(1), n) &= P^{S\_-}(B, n), \quad 0 \le n \le a-1, \end{aligned} \tag{32}$$

*and*

$$\begin{aligned} P^+(B, n) &= P^{S\_-}(B, n + b), \quad n \ge 1, \\ P^+(B, 0) &= \sum\_{n=a}^b P^{S\_-}(B, n). \end{aligned} \tag{33}$$

**Corollary 2.** *PS*− (*I*(*k*), *<sup>n</sup>*), 0 ≤ *n* < *a*, 1 ≤ *k* ≤ *c and PS*− (*<sup>B</sup>*, *<sup>n</sup>*), *n* ≥ 0 *satisfy the following equations:*

$$\begin{aligned} P^{S\_{-}}(I(k),n) &= \frac{P(I(k),n)}{1 - \sum\_{i=0}^{a-1} P(I(c),i)}, \quad 0 \le n \le a-1, 1 \le k \le c-1, \\\ P^{S\_{-}}(B,n) &= \frac{P(B,n)}{1 - \sum\_{i=0}^{a-1} P(I(c),i)}, \quad n \ge 0. \end{aligned} \tag{34}$$

**Proof.** When the service time distribution is exponential, service completions, real or potential, occur at random epochs. The probabilities, *PS*− (*I*(*k*), *<sup>n</sup>*), 0 ≤ *n* < *a*, 1 ≤ *k* ≤ *c* and *PS*− (*<sup>B</sup>*, *<sup>n</sup>*), *n* ≥ 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).
