2.3 General Behavior of Queueing Systems

In this section, we develop a simulation environment for queueing systems to explore their behavior and visualize queue-length dynamics.

The first aspect to address is system stability. Intuitively, when work arrives faster than it can be served, on average, the queue grows; if work arrives slower, the queue remains close to the origin. To illustrate this, Fig. 2.1 shows queue-length sample paths for several arrival rates and service capacities.

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_style("whitegrid")
tex_fonts = {
    "text.usetex": True,
    "font.family": "fourier",
    "axes.labelsize": 10,
    "font.size": 10,
    "legend.fontsize": 8,
    "xtick.labelsize": 8,
    "ytick.labelsize": 8,
}

plt.rcParams.update(tex_fonts)

Let \(a_k\) be the number of arrivals and \(c_k\) the service capacity in period \(k\). We write

We use Alg. 2.2.1 to compute the sample paths; jobs arriving in period \(k\) can be served in that same period. We simulate \(n=50\) periods with initial queue length \(L_0=40\) and mean arrival rate \(\lambda=6\). For each of ten independent sample paths, the arrivals \(a_k\) are uniformly distributed on \(\{5,6,7\}\), while the service capacity is constant at \(c_k=7\), and hence \(\mu=7\). The resulting stable queue-length trajectories are shown in the left panel of Fig. 2.1. For the right panel, we repeat the experiment with constant service capacity \(c_k=5\), and hence \(\mu=5\). Since \(\mu<\lambda\), the queue has positive drift and is unstable.

How about reducing the difference between \(\lambda\) and \(\mu\)? When \(\mu\) barely exceeds \(\lambda\), large fluctuations take longer to die out. To see this in quantitative terms, we run two simulations with smaller values for \(\mu\) while keeping \(\lambda=6\), and we compute the mean and variance of the number of jobs in the system.

For the first experiment, let \(c_k=6+B_k\), where \(B_k\sim\Bern{1/2}\) independently. Then \(\E{c_k}=6.5\). With \(n=5000\), \(L_0=0\), and the same arrival distribution as above, the simulation gives a sample mean queue length of \(0.45\) and a sample variance of \(0.64\).

For the second experiment, let \(c_k=6+B_k\), where \(B_k\sim\Bern{1/10}\) independently, so \(\mu=\E{c_k}=6.1\). With the same run length, initial queue length, and arrival distribution as above, the sample mean queue length is \(4.36\) and the sample variance is \(25.30\). Thus, reducing the excess service capacity from \(0.5\) to \(0.1\) substantially increases both the mean and the variance of the queue length.

Would variability affect the average queue length? We may suspect this because when \(a_{k}=1\) and \(c_{k}=1\) for all \(k\) and \(L_{0} = 0\), \(L_{k} = 0\) for all \(k\). If instead \(c_{i}=5\) when \(i\) is a multiple of \(5\) and \(0\) otherwise, queues do appear and disappear.¹¹Observe that the mean service time is still 1; only the periods in which service is offered changed. Let us use simulation to see how variability can influence the queue-length process. The goal is now to make Fig. 2.2.

We first need a way to quantify relative variability. Suppose the standard deviation of the service time is \(1\) minute. If the mean service time is 1 hour, we are inclined to call the service times very regular, while if the mean is just 1 minute, we would call it irregular. To capture this relative effect, we define the square coefficient of variation (scv) of the rv \(X\) as²²The scv is unit-less.

Continuing with the simulation, we seek a probability distribution for \(\{c_{k}\}\) with mean \(\mu\) and scv \(c^{2}\). A simple way to achieve this is to use a rv \(X\) such that \(\P{X=b} = p\) and \(\P{X=0} = 1-p\), and require that \(\E X = p b = \mu\) and \(C^2_X = c^{2}\). Using \(\V X = \E{X^{2}} - (\E X)^{2}\), we see that \(C^{2}_{X} = \E{X^{2}}/(\E X)^{2} - 1\). Rewriting, we get \(\E{X^{2}} = \mu^{2}(c^2+1)\), but we also know that \(\E{X^{2}} = p b^{2}\) for a Bernoulli-distributed rv. Solving for \(b\) and \(p\) in terms of \(\mu\) and \(c^{2}\) gives \(p = 1 / (1 + c^{2})\) and \(b=\mu/p = \mu (1+c^{2})\).

Thus, independent copies of \(X=bB\), with \(B\sim\Bern{p}\), provide a sequence with mean \(\mu\) and the required scv.

We simulate \(n=5000\) periods with \(L_0=0\), \(\lambda=6\), and \(\mu=6.5\), and compare scv values \(0.5\), \(1\), and \(2\). In the left panel of Fig. 2.2, arrivals are constant, \(a_k=6\), while the service capacities follow the two-point distributions constructed above. In the right panel, both arrivals and service capacities follow such distributions with the same scv. Applying the queue-length algorithm in each case shows that greater relative variability produces larger queue-length fluctuations.

In conclusion, tight capacity or high relative variability leads to large queues. To improve performance, we have basically four options: increase capacity, reduce load (or demand), reduce the scv of the arrivals or of the services.³³Of course, these options are not exclusive; for example, we can try to reduce the variability of the service process and add some capacity at the same time.