In this section, by using the renewal-reward theorem, we derive for the \(M^{B}/M/1\) queue expressions for \(\rho\), \(\E \W\), and \(\E \L\). To compute more difficult performance measures, such as the loss probability \(\P{\L>n}\), we need expressions for the stationary distribution \(\pi(n)\). We use level-crossing to derive a recursive scheme which, when combined with a policy to reject batches, can be used to compute these probabilities.

Assume that batches arrive as a Poisson process with rate \(\lambda_{B}\), and the number of jobs in a batch are iid \(\sim B\). Let \(S_{k,i}\) denote the service time of the \(i\)th job of batch \(k\); assume that the \(S_{k, i}\) are iid \(\sim S\) and independent of \(B\). By Eq. (4.1.2) the utilization \(\rho = \lambda_{B} \E B \E S\); as always, we require \(\rho<1\).

An arriving batch joins the end of the queue (if present). Once it arrives at the head of the queue, the batch moves in its entirety to the server¹¹ This is the same batching model as in Section 6.2. , and then the server processes the jobs one after another until the batch is finished. Writing \(\QQ^B\) for the number of batches in queue and \(\Ls^B\) for the number of jobs (if any) at the server, the average time \(\E{\W^{B}}\) a batch spends in queue is²² First the jobs at the server have to be served, then the ones in queue.

Substituting the equality³³ Little's law. \(\E{\QQ^B} = \lambda_{B} \E{\W^B}\) and simplifying yields

We can relate the average time \(\E \W\) a job spends in queue to \(\E{\W^{B}}\) by observing that \(\E{\Ls^B} + \E{\QQ^B} \E B = \E \L\),⁴⁴ The total number of jobs in the system is the number at the server plus the number in queue. The latter equals the number of batches in queue times the average batch size. so from Eq. (5.3.1)

where 1 follows from the memorylessness of the service time of the job in service (if any). Thus, to compute \(\E \W\) it remains to find an expression for \(\E{\Ls^B}\) in Eq. (5.3.2).

For this, we use the renewal reward theorem. Define⁵⁵ Much like the proof of Little's law. \(R(t) := \int_0^t \Ls^B(s) \d s\), where \(\Ls^{B}(s)\) is the number of jobs (of the batch in service) at the server at time \(s\), and \(Z_k = R(D_k)-R(D_{k-1})\) where \(D_k\) is the departure time of the \(k\)th batch. Suppose the \(k\)th batch has batch size \(B_{k}\), then

because \(B_{k}\) jobs are at the server during the service of the first job, \(B_{k}-1\) during the service of the second, \ldots Using that the \(S_{k,i} \sim S\),

As \(B_{k}\sim B\),

We now have all elements to fill in \(R=\alpha Z\). First, \(R = \lim_{t\to\infty}R(t)/t = \E{L_{S}^{B}}\). Second, we sample at batch departure moments, so \(\alpha=\delta_{B}\) which equals \(\lambda_{B}\) by rate stability. Third, \(Z=\E{Z_{k}}\). Combining all this, and using that \(\rho=\lambda_{B} \E B\E S\), we obtain

We can brush this up by realizing that⁶⁶ Recall that the scv of a rv \(B\) is \(C_B^2 = \V B/ (\E B)^{2}\). We assume that \(\E{B^{2}} < \infty\).

Substituting this into the above gives the final result

This formula is an important result. It tells us that \(\E \W\) increases when the load \(\rho\) or the variability \(C_{B}^{2}\) or \(\E B\) or \(\E S\) increases (even when the utilization \(\rho\) remains the same). As such, it gives us the most important clues about what to change if \(\E \W\)⁷⁷ Hence, by Little's law, \(\E \L\). has to be reduced.

Rather than using the time an entire batch spends at the server, we can concentrate on the expected time spent by a single job at the server. This will provide us with a simple expression for \(\P{\Ls^{B}=i}\), from which we also can derive Eq. (5.3.3).

In this case, let \(R_i(t) = \int_0^t \1{\Ls^{B}(s)=i} \d s\) be the total time there are \(i\) jobs at the server, so that \(R_i(t)/t \to \P{\Ls^{B}=i}\) as \(t\to\infty\). By sampling \(R(\cdot)\) at batch departure times \(\{D_k\}\), see Fig. 1, we see that

because only when \(B_{k}\geq i\), we can see \(i\) jobs at the server for the \(k\)th batch, and then \(S_{k,i}\) is the time this job occupies the server. With this,

since \(S_{k, i}\) and \(B_{k}\) are independent, and⁸⁸ Recall: \(\P{A} := \E{\1{A}}\). \(G(i-1) = \E{\1{B\geq i}} = \P{B> i-1}\) is the survivor function of \(B\). With the renewal reward theorem and rate stability, \(\delta_{B} = \lambda_{B}\), we conclude that

This expression provides us with another derivation of Eq. (5.3.4). Just note that \(\E{\Ls^{B}} = \sum_{i=0}^\infty i \P{\Ls^{B}=i} = \lambda_{B} \E S \sum_{i=1}^\infty i G(i-1)\) and use Opt 5.3.12 to simplify.

The second task is to use level-crossing to find a recursive formula for the state probabilities \(\{\pi(n)\}\).⁹⁹ Note: \(n\) means the number of jobs in the system, not the number of batches. For this we need to generalize our earlier level-crossing equation \(\lambda p(n) = \mu p(n-1)\) because when a batch arrives, multiple levels can be up-crossed, see Fig. 2. We use the PASTA property to express level-crossing in terms of \(\{\pi(n)\}\) instead of \(\{p(n)\}\).

The down-crossing rate of level \(n\) is easy: since jobs are still served one-by-one, we have \(\mu \pi(n+1)\).

For the up-crossings of level \(n\) we should consider batch arrivals that see \(m\leq n\) jobs in the system. Since batches arrive at a rate \(\lambda_{B}\), PASTA implies that \(\lambda_{B} \pi(m)\) is the arrival rate of batches that see \(m\) upon arrival.¹⁰¹⁰ That is, the total arrival process is thinned with a fraction \(\pi(m)\). When an arriving batch sees \(m\) jobs in the system, and level \(n\) has to be up-crossed, the batch must contain more than \(n-m\) jobs. As \(G(n-m)\) is the probability that the batch size is larger than \(n-m\) and as batch sizes are independent of arrival times, such batches arrive as a Poisson process with rate \(\lambda_{B} \pi(m) G(n-m)\). Finally, to compute the rate at which level \(n\) is up-crossed, we merge all these thinned Poisson processes into a new Poisson process with rate \(\lambda_{B} \sum_{m=0}^n \pi(m) G(n-m)\).

By level-crossing, the up-crossing and down-crossing rates must match, thus,

It is an important check¹¹¹¹ Opt 5.3.13. to retrieve Eq. (5.3.4) from the expression \(\E\L = \sum_{n} n \pi(n)\) and Little's law.

It remains to find the normalization constant for \(\pi\). However, this is not possible in general because Eq. (5.3.5) does not lead to a closed form expression for \(\pi(n)\).¹²¹² Compare Eq. (5.1.1). Numerically, though, we can start with \(\pi(0)=1\), apply the recursion, and truncate when the (unnormalized) probabilities become negligible. The RV class then normalizes the result. The function mxm1 applies this idea.

from random_variable import NumericRV as RV

labda_B, mu = 1, 3
B = RV({1: 1, 2: 1, 3: 1})
rho = labda_B * B.mean() / mu
rho_B = labda_B / mu


def mxm1(eps=1e-10):
    "Compute stationary distribution of number of jobs in the system."
    if rho >= 1:
        print("The load is too high")
        quit()
    pi, n = [1], 0
    batch_sizes = [int(m) - 1 for m in B.support()]
    while n < batch_sizes[-1]:
        res = sum(pi[m] * B.sf(n - m) for m in range(n + 1))
        pi.append(res * rho_B)
        n += 1
    while pi[-1] > eps:
        res = sum(pi[n - m] * B.sf(m) for m in batch_sizes)
        pi.append(res * rho_B)
        n += 1
    return RV({i: pi[i] for i in range(len(pi))})

A different approach is to completely bypass the normalization problem by blocking the demand above some level \(K\). There are three common policies for deciding which jobs in a batch to accept.

1. Complete acceptance: accept all batches that arrive when the system contains \(K-1\) or fewer jobs, and reject the entire batch when the system contains \(K\) or more jobs. For this case, we make the assumption that the batch size is bounded by the number \(K_{B}\) so that \(L(t) \leq K + K_{B}\).¹³¹³ Ex 5.3.8.
2. Partial acceptance: accept whatever fits of a batch, i.e., up to \(K\), and reject the rest.¹⁴¹⁴ Ex 5.3.9.
3. Complete rejection: if a batch does not fit entirely into the system, reject it completely.¹⁵¹⁵ Ex 5.3.10.

In Ex 5.3.9--Ex 5.3.10 we will derive recursions similar to Eq. (5.3.5) to compute \(\{\pi(n)\}\) for these different blocking models.

def complete_rejection(K):
    pi, n = {}, 0
    pi[0] = 1
    while n < K:
        pi[n + 1] = sum(
            pi[m] * (S.sf(n - m) - S.sf(K - m)) for m in range(n)
        )
        pi[n + 1] *= rho
        n += 1
    return RV(pi)

Suppose \(B\) is geometrically distributed¹⁶¹⁶ Sometimes this distribution is called the First Success distribution. There is no consensus on the definition of the geometric distribution, see Wikipedia. with \(\E B = 1/p\).

def complete_acceptance(K):
    pi, n = {}, 0
    pi[0] = 1
    while B.sf(n - K) > 0:
        pi[n + 1] = sum(pi[m] * B.sf(n - m) for m in range(min(n, K) + 1))
        pi[n + 1] *= rho_B
        n += 1
    return RV(pi)

def partial_acceptance(K):
    pi, n = {}, 0
    pi[0] = 1
    while n < K:
        pi[n + 1] = sum(pi[m] * B.sf(n - m) for m in range(n + 1))
        pi[n + 1] *= rho_B
        n += 1
    return RV(pi)

def complete_rejection(K):
    pi, n = {}, 0
    pi[0] = 1
    while n < K:
        pi[n + 1] = sum(
            pi[m] * (B.sf(n - m) - B.sf(K - m)) for m in range(n + 1)
        )
        pi[n + 1] *= rho_B
        n += 1
    return RV(pi)