We start with a case to provide motivation for studying queueing systems. Then we develop a set of fundamental recursions to simulate the case and evaluate the efficacy of several policies to improve the system.

At a mental-health department, five psychiatrists perform intakes¹¹ Interviews with new patients. to determine the best treatment process once patients enter the system. There are complaints about the time patients have to wait for the intake; the desired waiting time is around two weeks, but the realized waiting time is sometimes more than three months. This is unacceptably long, but what to do about it?

The five psychiatrists put forward various suggestions to reduce waiting times.

1. Give all psychiatrists the same weekly capacity to perform intakes. Motivation: Currently, one psychiatrist performs 9 intakes per week, one performs 3, and the other three psychiatrists only perform 1. This is not a problem during weeks when all psychiatrists are present; however, psychiatrists take vacations, visit conferences, etc. Hence, if the psychiatrist with the highest weekly number of intakes goes on leave, the total intake capacity is significantly reduced.
2. Synchronize holidays.²² Later, with the models of Chapter 6, we will see that this will not work, in fact, synchronization might be detrimental. Motivation: to reduce the variation in the service capacity, the psychiatrists plan their holidays consecutively. However, perhaps it is better to work at full capacity or not at all.

After consultation, another suggestion was proposed:

3. Control the intake capacity as a function of the waiting time.

Motivation: in analogy with supermarkets, increase (decrease) capacity when the queue is long (short). Observe that people often object to such policies because they believe that they have to do more work. However, this reasoning is incorrect when the arrival rate does not depend on the queue length. To see this, suppose that 1000 patients per year need treatment. This number does not depend on whether patients spend time in the queue or not.

We next develop a method to simulate the behavior of the system over time so that we can evaluate the effect of the above suggestions. We use simulation because mathematical analysis is intractable.³³ In case you doubt this, try to analyze transient multi-server queueing systems with vacations, of which this is an example.

Let us start with analyzing the two suggestions of the psychiatrists by making a queueing model. Often, the easiest way to construct queueing processes is to divide time into periods of equal length, and develop recursions that capture the behavior of the queue from one period to the next.⁴⁴ Using fixed-length periods has the advantage that we do not have to specify inter-arrival times or service times of individual customers. In this case, only the number of arrivals in a period and the number of potential services are relevant.

Let us define:

• \(a_k =\) the number of jobs that arrive in period \(k\),
• \(c_k =\) the capacity, i.e., the maximal number of jobs that can be served, during period \(k\),
• \(d_k =\) the number of jobs that depart just before the end of period \(k\),
• \(\L_k =\) the system size, i.e., the number of jobs in the system at the end of period \(k\).

In what follows, we also call \(a_k\) the size of the batch arriving in period \(k\). Note that the definition of \(a_k\) is subtle: we may assume that the arriving jobs arrive either at the start or at the end of the period. In the first case, the jobs can be served in period \(k\); in the latter case, they cannot be served in period \(k\).

Since \(\L_{k-1}\) is the length of the system at the end of period \(k-1\), it must also be the number of customers at the start of period \(k\). Assuming that jobs arriving in period \(k\) cannot be served in period \(k\), the number of customers \(d_{k}\) that depart in period \(k\) and the number that remain in the system \(\L_{k}\) at the end of period \(k\) are given by the recursions

Thus, if we are given \(\L_0\) and arrivals \(\{a_{k}\}_{k=1}^{n}\) and capacities \(\{c_{k}\}_{k=1}^{n}\), we can obtain \(\{d_{k}\}_{k=1}^{n}\) and \(\{L_{k}\}_{k=0}^{n}\).

Typically, we can use company data to model the arrival process \(\{a_k\}\) and the capacity \(\{c_k\}\), and feed these data into a computer to perform recursions Eq. (2.1.1). If insufficient data is available, we build a probability model for \(\{a_{k}\}\) and \(\{c_{k}\}\) and generate random deviates with these distributions.

Another interesting recursion is this:⁵⁵ \(x^{+} := \max\{x, 0\}\), \(x^{-} := [-x]^{+}\).

Ex 2.1.9 asks you to prove that this generates a system in which jobs can be served⁶⁶ If there is sufficient capacity. in the period they arrive.

Continuing with the case of the five psychiatrists, we can now analyze how the first and second suggestion would affect the waiting times.

As a first step in the analysis, we model the arrival process of patients as a Poisson process. The duration of a period is taken to be a week. The average number of arrivals per period, based on company data, is slightly less than \(12\) per week; in the simulation we set \(\lambda= 11.8\) per week.

We model the capacity in the form of a matrix such that row \(i\) corresponds to the weekly capacity of psychiatrist \(i\):

Thus, psychiatrists 1, 2, and 3 perform one intake per week, the fourth performs 3, and the fifth performs 9. The sum over column \(k\) is the total service capacity for week \(k\) of all psychiatrists together.

With such a matrix, it is easy to construct other capacity schemes. A more balanced scheme would be like this:

We include the effects of holidays by setting the capacity of one or some psychiatrists to \(0\) in a certain week. Let us assume that every week just one psychiatrist is on leave such that their holiday schemes rotate. To achieve this, we set the entries of the matrix like this:

Clearly, under this scheme, the average capacity is \(12\) patients per week. There is another simple holiday scheme: when all psychiatrists take holiday in the same week:

Of course, random holiday schemes are more likely, but these two are perhaps the simplest that result in the same average capacity.

By also applying the same holiday plans to the unbalanced case of Eq. (2.1.3), we obtain four different capacity schemes: unbalanced versus balanced capacity, and spread out holidays versus simultaneous holidays. We can now use the recursions Eq. (2.1.1) to simulate the queue length process.

The results are shown in Fig. 1. We plot, for each period, the largest and smallest queue that occurred under all four capacity plans that result from following the first and second suggestions of the psychiatrists. The graphs show that these suggestions hardly affect the behavior of the queue-length process.

This desk study brings us to an important observation: Suggestions 1 and 2 are unlikely to have a positive impact on waiting times. So, rather than setting up an experiment in real life, we should look for other rules.

Now we consider Suggestion 3, which consists of doing more intakes when it is busy and fewer when it is quiet. A simple rule is to let the capacity for week \(k\) depend on the queue length at the end of the previous week, for instance,

We can take \(e=1\) or \(2\), or perhaps a larger number; the larger \(e\), the larger the control we exercise.⁷⁷ Why would it not make much sense to set the lower threshold, here 24, smaller than 12? Hint, do we use the capacity of 12 patients efficiently? We can, of course, also adapt the thresholds 48 and 24.

Consider three different control levels, \(e=1\), \(e=2\), and \(e=5\); when \(e=1\), one psychiatrist does one additional intake, when \(e=5\), all do an extra intake in that week. The right panel of Fig. 1 shows a striking result. The queue does not explode anymore; even taking \(e=1\) has a large effect.

In conclusion, this simulation experiment shows that changing holiday plans or spreading the work more evenly over multiple servers, i.e., psychiatrists, may not significantly affect queueing behavior. However, controlling the service rate as a function of the queue length can dramatically improve the situation.

We next discuss some general points that need further attention. In Eq. (2.1.1) we assume that the jobs that arrive in period \(k\) cannot be served in period \(k\). If the situation is such that jobs can be served in the same period as they arrive, then Eq. (2.1.1) should be changed to⁸⁸ Ex 2.1.9.

Which of Eq. (2.1.1) or Eq. (2.1.8) to choose depends on what we need to model; in general, no rule is `best', and what is `good' depends essentially on the context. For instance, if we want to be `on the safe side', then it is perhaps best to use Eq. (2.1.1) because with this rule, we overestimate the queue lengths, while with Eq. (2.1.8) we underestimate the queue lengths.

Note next that in the computation of \(d_{k}\) we make a fundamental modeling choice: If there are jobs in the system and the server capacity \(c_{k} > 0\), the server will serve jobs. However, there is no formal need to serve jobs, even when there is service capacity available. A simple reason not to serve jobs is when it is too expensive; for example, if there is very little demand for a flight on some day, then the flight can be canceled. In the models we consider, we will not allow jobs to wait when there is capacity available; in other words, our service processes are assumed to be work-conserving.

The above recursions only construct \(\{d_{k}\}\) and \(\{L_k\}\). If we also need information about the sojourn times, that is, the time jobs spend in the system, it is necessary to specify the service discipline, i.e., a scheduling rule that decides on the order in which jobs in queue are taken into service. In this book, we henceforth assume that jobs move to the server in the order in which they arrive. This is known as First-In-First-Out (FIFO).⁹⁹ First-Come, First-Served (FCFS) is also often used. There are many other scheduling rules, such as Last-In-First-Out; we do not discuss these here.

It is quite remarkable that with Eq. (2.1.1) the computation of the system length process \(\{L_{k}\}\) is very simple, but it is much harder to compute the sojourn time \(\J_{k}\) for the jobs arriving in period \(k\). In fact, the best we can obtain are bounds on the sojourn time \(\J_{k}\) such that \(\J_k^{-} \leq \J_{k}\leq \J_k^{+}\), where

To see how this works, suppose that \(\L_{k-1}=20\), \(a_{k}=6\), and \(c_k= c= 2\) for all \(k\). Any job that arrives in the \(k\)th period, must wait at least \(L_{k-1}/c = 20/2 = 10\) periods just to get access to the server under the FIFO discipline, and we can be sure that the last of these jobs is served after \((20+6)/2=13\) periods.

We note that there is a difference between waiting time \(\W\) and sojourn time \(\J\). The former is the time spent by jobs in the queue, the latter is the time in the system, which is the waiting time plus the time at the server. However, in Eq. (2.1.1) and Eq. (2.1.8) jobs are not really in service; we only count the jobs in the system at the end of a period. Thus, in this model, the number of jobs in the system coincides with the number in queue.

Finally, once the simulation is finished, we compute some performance measures to analyze how the (simulation of the) queueing system behaved. Besides making graphs, we can be interested in the arrival rate, the capacity¹⁰¹⁰ Also known as the service rate. and mean number in the system,

Similarly, we can take the averages over \(d_{k}\), \(\J_k^{-}\) and \(\J_k^{+}\). The variances of these (simulated) rvs can also be of interest. Perhaps the most important measure to check is the load, which is the rate at which work arrives per server. For instance, if \(\lambda=5\) jobs arrive per hour, and each job requires 20 minutes of work (on average), then we say that the load is \(5\times 20/60 = 5/3\). Since one server can do only 1 hour of work per hour, we need at least two servers to deal with this number of jobs. For a simulation in discrete time, we define the load as

We expect that when \(\rho>1\), long queues build up.¹¹¹¹ In Section 2.3 we discuss the concept of load more formally.

In general, with recursions like Eq. (2.1.1) we can carry out numerous what-if analyses.

For example, a hospital considers buying a second MRI scanner because the current one is saturated: it is used from 8 am to 6 pm.¹²¹² That is, 10 h per working day. If the hospital can identify some additional capacity, it might postpone the purchase of an additional MRI scanner for a few years. Suppose that a percentage of the staff is willing to work from 8 pm to 11 pm, say, and that 30\% of the patients¹³¹³ Perhaps some staff and patients even prefer to work/have scans in the evening. are ready to come to the hospital for a scan in the evening, then the capacity would increase by about \(3/10 = 30\%\). To see whether this idea is interesting, we can use Eq. (2.1.1) to let the computer make a graph of the influence on \(L\).

Recursions are not just useful for modeling queueing systems. For instance, a common model in population dynamics has this form

where \(N_t\) is the population size at year \(t\), \(B_{t}\) the number of births, \(I_{t}\) the immigration, \(D_{t}\) the deaths, and \(E_{t}\) the emigration.

Yet another use of recursion is to understand how the reflection of infrared light by carbon dioxide affects the temperature of the earth. The atmosphere is modeled as a sequence of layers so that, per layer, the temperature, air density, and chemical composition are considered constant. These properties determine how much infrared light, which is initially reflected by the surface of the earth, reaches outer space. In each of the layers, we need to specify what fraction of infrared photons are transmitted to a higher layer, absorbed, or reflected again to a lower level. Thus, this physical process can be modeled as a queueing system in which photons (customers) go from one layer (station) to another layer (station) or are absorbed (leave the network).¹⁴¹⁴ https://github.com/scienceetonnante/RadiativeForcing/blob/main/RadiativeForcing.py.

Formulating recursions is the key step in making useful models; you should practice this a lot. In some exercises the solution is stated in code so that you practice reading (and thinking in) code.

d1[k] = min(L1[k - 1], c[k])
L1[k] = L1[k - 1] - d1[k] + a1[k]
c2[k] = min(L2[k - 1], c[k] - d1[k])
d2[k] = min(L2[k - 1], c2[k])
L2[k] = L2[k - 1] - d2[k] + a2[k]

c2[k] = min(L2[k - 1], r2)
d1[k] = min(L1[k - 1], c[k] - c2[k])
L1[k] = L1[k - 1] - d1[k] + a1[k]
d2[k] = min(L2[k - 1], c[k] - d1[k])
L2[k] = L2[k - 1] - d2[k] + a2[k]

The next true-false question relates to the control rule Eq. (2.1.7).

if L[k - 1] >= 48:
    c[k] = 12 + e
elif L[k - 1] <= 24:
    c[k] = 12 - e

For real-life situations the recursions can become quite intricate.¹⁵¹⁵ This example is a somewhat simplified case from a bachelor thesis.

Ex 2.1.7

A nurse takes blood samples in two departments of a hospital. It takes time to walk from one location to another. The nurse serves all patients in one location, then walks to the other location, serves all patients there, walks back to the first location, and so on.

Give the recursions to see how this switching behavior affects the queue lengths at both locations.

Hint

The recursions can be found by introducing an extra variable \(p_k\) that specifies the production state of the nurse. If \(p_k=0\), the nurse serves the first queue, if \(p_k=1\) the nurse moves from queue 1 to queue 2, if \(p_k=2\) the nurse serves queue 2. Finally, if \(p_k=3\) the nurse moves from queue 2 to queue 1. Thus, the production state cycles from \(0\to 1\to 2\to 3 \to 0\).

Solution

Did you actually try? Maybe see the 'hints' above!:

Solution, for real

Recall that code like (3 <= 4) implements the indicator function. Focus on the logic of for loop, and on how the state of the nurse is encoded. The rest is just code to make a nice figure. BTW, making nice visual reports is an important skill (but not so much for this course).

Python

Copy

import numpy as np
import matplotlib.pyplot as plt
from latex_figures import apply_figure_settings
from enum import Enum, auto

apply_figure_settings(use=True)


class State(Enum):
    A_TO_B = auto()
    B_TO_A = auto()
    AT_A = auto()
    AT_B = auto()


rng = np.random.default_rng(3)
thres_A, thres_B = 40, 10

num_weeks = 100
a1 = rng.poisson(3, size=num_weeks)
a2 = rng.poisson(4, size=num_weeks)
c1 = rng.poisson(7, size=num_weeks)
c2 = rng.poisson(8, size=num_weeks)

LA = np.zeros(num_weeks)
LB = np.zeros(num_weeks)
LA[0], LB[0] = 30, 50

p = np.zeros(num_weeks, dtype='object')  # production state
p[0] = State.AT_A

for k in range(1, num_weeks):
    if p[k - 1] == State.AT_A and LA[k - 1] <= thres_A:
        p[k] = State.A_TO_B
    elif p[k - 1] == State.A_TO_B:
        p[k] = State.AT_B
    elif p[k - 1] == State.AT_B and LB[k - 1] <= thres_B:
        p[k] = State.B_TO_A
    elif p[k - 1] == State.B_TO_A:
        p[k] = State.AT_A
    else:  # don't change
        p[k] = p[k - 1]

    d1 = min(LA[k - 1] + a1[k], c1[k] * (p[k] == State.AT_A))
    d2 = min(LB[k - 1] + a2[k], c2[k] * (p[k] == State.AT_B))
    LA[k] = LA[k - 1] + a1[k] - d1
    LB[k] = LB[k - 1] + a2[k] - d2


# Convert cm → inches for Matplotlib
cm_to_inch = 1 / 2.54
fig_width = 6 * cm_to_inch
fig_height = 5.5 * cm_to_inch

xx = range(num_weeks)
plt.figure(figsize=(fig_width, fig_height))
plt.step(xx, LA, ls=":", lw=0.7, where="pre", label="LA")
plt.step(xx, LB, ls=":", lw=0.7, where="pre", label="LB")
plt.xlabel('Time')
plt.legend()
plt.tight_layout()
plt.savefig("/figures/nurse.png", dpi=300)

import numpy as np
from numpy.random import default_rng

rng = default_rng(3)
num_weeks = 1000
a1a = rng.poisson(10, num_weeks) # new arriving cases
a2 = np.zeros(len(a1a))          # arrivals at judge
a1b = np.zeros(len(a1a))         # arriving for a second step at the clerks
d1a = np.zeros(len(a1a))         # departures from first round at clerks
d2 = np.zeros(len(a1a))          # departures from judge
d1b = np.zeros(len(a1a))         # departures after second round at clerks
L1a = np.zeros(len(a1a))         # amount of work in hours at the first stage at clerks
L2 = np.zeros(len(a1a))          # amount of work in hours at judge
L1b = np.zeros(len(a1a))
c1 = 80 * np.ones(len(a1a))      # weekly hours available of clerks
c2 = 34 * np.ones(len(a1a))      # hours available of judges

ES1a = 5                         # work in hours per case at first stage
ES1b = 2                         # work in hours per case at second stage
ES2 = 3                          # work of judge per case


for k in range(num_weeks):
    d1a[k] = min(L1a[k - 1], c1[k])
    L1a[k] = L1a[k - 1] + a1a[k] * ES1a - d1a[k]
    a2[k] = d1a[k] // ES1a
    d2[k] = min(L2[k - 1], c2[k])
    L2[k] = L2[k - 1] + a2[k] * ES2 - d2[k]
    a1b[k] = d2[k] // ES2
    d1b[k] = min(L1b[k - 1], c1[k] - d1a[k])
    L1b[k] = L1b[k - 1] + a1a[k] * ES1b - d1b[k]

print(f"{L1b.mean()=:.2f}")
print(f"{L2.mean()=:.2f}")
print(f"{L1a.mean()=:.2f}")