Consider a sample path of \(\{\L(t) : t\geq 0\}\) of a queueing or inventory system in which jobs arrive and depart one by one. We say that the sample path up-crosses level \(n\) at time \(t\) when the number of jobs in the system changes from \(n\) to \(n+1\) due to an arrival, that is, when \(L(t-)=n\) and \(L(t)=n+1\). The sample path down-crosses level \(n\) at time \(t\) when, due to a departure, \(L(t-)=n+1\) and \(L(t)=n\). Clearly, the number of up-crossings and down-crossings cannot differ by more than 1 at any time, because a down-crossing of level \(n\) can only occur after an up-crossing (or the other way around). This simple idea is a key tool in the analysis of stationary stochastic systems.

Similar to the definitions for \(A(n,t)\) and \(\lambda(n)\), let¹¹ Note that we write \(L(D_{k})\) and not \(L(D_{k}-)\): to leave \(n\) jobs behind, the system must contain \(n+1\) jobs just prior to the departure.

\begin{align*} D(n,t) &:= \sum_{k=1}^{D(t)} \1{\L(D_k) = n}, & \mu(n+1) := \lim_{t\to\infty} \frac{D(n,t)}{Y(n+1,t)}, \tag{4.6.1} \end{align*}

denote the number of departures up to time \(t\) that leave \(n\) jobs behind, and the departure rate from state \(n+1\).

Proof

Since jobs arrive one by one, \(|L(t) - L(t-)| \leq 1\) for all \(t\geq 0\). This implies that

\begin{equation*} |A(n,t) - D(n,t)| \leq 1, \tag{4.6.3} \end{equation*}

hence,

\begin{equation*} \lim_{t\to\infty} \frac{A(n,t)}t = \lim_{t\to\infty} \frac{D(n,t)}t. \tag{4.6.4} \end{equation*}

The claim follows by noting that

\begin{align*} \lim_{t\to\infty} \frac{A(n,t)}t &= \lim_{t\to\infty} \frac{A(n,t)}{Y(n,t)}\frac{Y(n,t)}t = \lambda(n) p(n),\\ \lim_{t\to\infty} \frac{D(n,t)}t &= \lim_{t\to\infty} \frac{D(n,t)}{Y(n+1,t)}\frac{Y(n+1,t)}t = \mu(n+1) p(n+1). \end{align*}

Suppose we can specify²² In Section 5.1 and onward we will model many queueing situations by making suitable choices for \(\lambda(n)\) and \(\mu(n)\). the arrival and service rates \(\lambda(n)\) and \(\mu(n)\), then we can easily compute the long-run fraction of time \(p(n)\) that the system contains \(n\) jobs. To see this, rewrite Eq. (4.6.2) as

\begin{equation*} p(n+1) = \frac{\lambda(n)}{\mu(n+1)}p(n). \tag{4.6.5} \end{equation*}

A simple recursion yields:

\begin{equation*} p(n+1) = \frac{\lambda(n)\lambda(n-1)\cdots \lambda(0)}{\mu(n+1)\mu(n)\cdots \mu(1)}p(0). \end{equation*}

Since \(p(n)\), \(n\geq 0\), represent probabilities,

\begin{equation*} 1 = \sum_{n=0}^\infty p(n)= p(0)\left(1+\sum_{n=0}^\infty \frac{\lambda(n)\lambda(n-1)\cdots\lambda(0)}{\mu(n+1)\mu(n)\cdots \mu(1)}\right). \end{equation*}

Thus, we can normalize \(p(n)\) by setting

\begin{equation*} p(0)^{-1} = 1+\sum_{n=0}^\infty \frac{\lambda(n)\lambda(n-1)\cdots\lambda(0)}{\mu(n+1)\mu(n)\cdots \mu(1)} \tag{4.6.6} \end{equation*}

Finally, once we have \(p(n)\), it is easy to compute the time-average number of jobs in the system and the long-run fraction of time it contains at least \(n\) jobs:

\begin{align*} \E\L &= \sum_{n=0}^\infty n p(n), & \P{\L \geq n} &= \sum_{i=n}^\infty p(i). \end{align*}

With the definitions of \(A(n,t)\) and \(D(n,t)\), we establish a relation between \(\pi(n)\) and \(\delta(n)\), i.e., statistics as obtained by arrivals and departures. Define, analogous to \(\pi(n)\),

\begin{equation*} \delta(n) := \lim_{t\to\infty} \frac{D(n,t)}{D(t)} \end{equation*}

as the long-run fraction of jobs that leave \(n\) jobs behind. From Eq. (4.6.4) and using that \(A(n,t) \approx D(n,t)\),

\begin{equation*} \frac{A(t)}t \frac{A(n,t)}{A(t)} = \frac{A(n,t)}t \approx \frac{D(n,t)}t = \frac{D(t)}t \frac{D(n,t)}{D(t)}. \end{equation*}

Taking limits at the left and right, and using Eq. (4.1.1), we obtain for the \(G/G/c\) queue³³ Because customers arrive and leave as single units in a (rate-stable) \(G/G/c\) queue.

\begin{equation*} \lambda \pi(n) = \delta \delta(n). \tag{4.6.7} \end{equation*}

Hence,

\begin{equation*} \lambda = \delta \iff \pi(n) = \delta(n). \tag{4.6.8} \end{equation*}

When jobs can be blocked, we must distinguish between jobs that arrive and those that actually enter.⁴⁴ Blocked jobs do not enter, hence do not produce up-crossings. To handle this, define

\begin{align*} A^{+}(n,t) &:= \sum_{k=1}^{A(t)} \1{L(A_k-) = n}\1{\L(A_k) \geq n+1}, & \lambda^{+}(n) &:= \lim_{t\to\infty} \frac{A^{+}(n,t)}{Y(n,t)}, \tag{4.6.9} \end{align*}

i.e., \(A^{+}(n,t)\) counts arrivals that up-cross level \(n\), and \(\lambda^{+}(n)\) corresponds to the rate at which jobs enter when the system contains \(n\) jobs.⁵⁵ \(\lambda(n) > \lambda^{+}(n)\) when a fraction of jobs is blocked, otherwise \(\lambda(n) = \lambda^{+}(n)\). We should replace \(\lambda(n)\) by \(\lambda^{+}(n)\) at all other relevant places, for instance, Eq. (4.6.2) becomes

\begin{equation*} \lambda^{+}(n) p(n) = \mu(n+1) p(n+1). \tag{4.6.10} \end{equation*}

It is important to realize that level-crossing arguments in the simple form as presented above cannot always be applied, since it is not always possible to split the state space into two disjoint subsets by `drawing a line' between two states. For a more general approach, we focus on a single state and count how often this state is entered and left.

Specifically, define \(I(n,t) = A(n-1,t) + D(n,t)\) as the number of times the queueing process enters state \(n\) either by an arrival from state \(n-1\) or a departure leaving \(n\) jobs behind. Similarly, \(O(n,t) = A(n,t) + D(n-1,t)\) counts how often state \(n\) is left either by an arrival (to state \(n+1\)) or a departure (to state \(n-1\)). Just like Eq. (4.6.3), it is evident that \(|I(n,t)-O(n,t)|\leq 1\), and this implies⁶⁶ Ex 4.6.10.

\begin{equation*} \lambda(n-1)p(n-1)+\mu(n+1)p(n+1) = (\lambda(n)+\mu(n))p(n). \tag{4.6.11} \end{equation*}

These equations hold for any \(n\geq 1\) and are known as the balance equations. These equations are useful for analyzing queueing networks.

Level-crossing arguments can also be applied to continuous systems. For instance, we can find a closed-form expression for the waiting time of the \(M/G/1\) queue. However, this is a bit more complicated as the argument hinges on an integral equation. If you do not like developing mathematical skills, you may skip this derivation.

Suppose that the waiting time \(W\) of the \(M/G/1\) queue has a density \(f\), that is, the cdf \(F(x) = \P{W\leq x}\) has a derivative \(f(x)\) for all \(x>0\).⁷⁷ A formal proof would first establish that \(f\) exists. ^{, 88} At \(x=0\), \(F\) is only differentiable from the right because \(\P{W=0}>0\), while \(\P{W\leq x} = 0\) for \(x<0\). Consider some level \(x>0\). By PASTA, the rate of arriving jobs that see a waiting time \(y \in [0, x)\) is equal to \(\lambda f(y)\).⁹⁹ Observe that we use the splitting property of Poisson processes. To up-cross level \(x\) from state \(y< x\), the service time must be at least \(x-y\), the probability of which is \(G(x-y) = \P{S>x-y}\). Adding up the rates for all possible waiting times below \(x\) gives that \(\lambda\int_0^x f(y) G(x-y)\d y\) is the up-crossing rate of level \(x\). The down-crossing rate equals \(f(x)\) because the server drains work at unit rate, so level \(x\) is crossed downward whenever the workload is at \(x\). Equating these rates gives the integral equation \(\lambda \int_0^{x} f(y) G(x-y)\d y = f(x)\). For \(x=0\), we need to consider the atom \(\P{W=0}\). In that case, \(\lambda \P{W=0} = f(0+)\). We can combine these observations into the following level-crossing equation:

\begin{equation*} \lambda\int_0^x G(x-y) \d F(y) = f(x). \tag{4.6.12} \end{equation*}

Using this integral equation

\begin{align*} \E W &= \int_0^{\infty} x f(x) \d x = \lambda \int_0^{\infty} x \int_0^{x} G(x-y) \d F(y) \d x \\ &= \lambda \int_0^{\infty} \int_0^{\infty} x \1{y \leq x} G(x-y) \d x \d F(y) \\ &= \lambda \int_0^{\infty} \int_0^{\infty} (u+y) G(u) \d u \d F(y) \\ &\stackrel 1 = \lambda \int_0^{\infty} (\E{S^2}/2 + y \E S) \d F(y) \\ &\stackrel 2= \lambda \E{S^2}/2 + \lambda \E S \E \W. \end{align*}

Step 1 follows from substituting \(u=x-y\), and 2 from the fact that \(\E \W = \int_{0}^{\infty} y \d F (y)\). Bringing the term \(\lambda \E S \E \W = \rho \E \W\) to the left, and dividing by \(1-\rho\) gives

\begin{align*} \E \W &= \lambda \frac{\E{S^2}}{2(1-\rho)}. \end{align*}

As a last point of interest, let us derive the density of the waiting time for the \(M/M/1\) queue. As in this case \(S\sim \Exp{\mu}\), \(f\) must satisfy

\begin{align*} f(x) = \lambda \int_0^x e^{-(x-y)\mu} \,\d F(y) = \lambda \int_0^x e^{-(x-y)\mu} f(y) \d y + \lambda e^{-x \mu} F(0). \end{align*}

Multiplying both sides by \(e^{\mu x}\) and defining \(g(x) = e^{\mu x} f(x)\), this can be rewritten to \(g(x) = \lambda \int_0^x g(y) \d y + \lambda F(0)\). Differentiating the LHS and RHS with respect to \(x\) gives the differential equation \(g'(x) = \lambda g(x)\). Hence, \(g(x) = A e^{\lambda x}\), and therefore, \(f(x) = A e^{-(\mu-\lambda)x}\). At \(x=0\), \(f(0) = \lambda \P{W=0} = \lambda p(0)\). Since \(p(0)=1-\rho\), it follows that

\begin{equation*} f(x) = \lambda(1-\rho) e^{- \mu(1-\rho)x}, \quad x > 0. \tag{4.6.13} \end{equation*}

Another way to obtain this result is to condition on the number \(\L\) of jobs found upon arrival, and then using that the total service of these \(\L\) customers has a gamma density:¹⁰¹⁰ If \(L(t) = k > 0\), the new arrival finds one job at the server and \(k-1\) in queue. Each of these jobs in the system has to be served before it is the turn of the new job.

\begin{align*} f(x) &= \sum_{k=1}^\infty \P{W = x | \L = k}\P{\L = k } \\ &= \sum_{k=1}^\infty \P{S_1 + S_2 + \cdots + S_k = x}\P{\L = k } = (1-\rho) \sum_{k=1}^\infty \rho^k \mu \frac{(\mu x)^{k-1}}{(k-1)!}e^{-\mu x}\\ &= (1-\rho) \lambda e^{-\mu x}\sum_{k=0}^\infty \rho^{k} \frac{(\mu x)^{k}}{k!} = (1-\rho) \lambda e^{-\mu (1-\rho)x}. \end{align*}