Case Study: Multiple-Server Queue - Modelling and Simulation Cookbook

We model the behaviour of the queue of customers waiting at the cashier desks in a supermarket.

Topics:

Discrete Event Simulation

Objectives:

How does the practical implementation of a discrete-event model work from scratch?
What are performance bottlenecks of the implementation and how can they be improved?
What is the impact of randomess on the results?

import datetime
from abc import abstractmethod
from collections.abc import Callable
from heapq import heapify, heappop, heappush
from typing import Self

import numpy as np
%matplotlib ipympl
from matplotlib import pyplot as plt

Motivation / System¶

A modeller needs to model the behaviour of a single supermarket queue with $k$ cashiers. He finds that a new customers arrives, on average, every $\mu_a$ seconds and that it takes on cashier, on average, $\mu_s$ seconds to serve a customer.

Conceptual Model¶

The modeller opts for a classic discrete event model and specifies the model using the Event Graph below. The times $t_a$ and $t_s$ are hereby random variables with $E(t_a)=\mu_a$ and $E(t_s)=\mu_s$ .

Implemented Model¶

Below, we define a base-class for a generic discrete event model from scratch. The class is driven by event notices that extend the Event class. These object need to implememt the eval which executes what happens when they are executed.

# from heapq import heapify, heappush, heappop  # improves performance by using a heap for the event queue
class States:
    pass


class Event:
    """Base class for an event-notice"""

    @abstractmethod
    def eval(self, states: States) -> list[tuple[Self, float]]:
        """What happens when the event takes place.
        :param states: some object carrying the current value of the states
        :return: list of tuples with new Events and corresponding delays to be scheduled by the simulation
        """
        return list()


class DESSimulator:
    def __init__(self):
        """Base implementation of a discrete event simulator"""
        self.t = 0  # internal time
        self.event_list: list[tuple[float, Event]] = list()  # internal event list
        # heapify(self.event_list)

    def add(self, event: Event, delay: float):
        """Schedules a new event to the list with a given delay in the future
        :param event: event notice
        :param delay: delay greater or equal to 0
        :return:
        """
        if delay < 0:
            return
        self.event_list.append((self.t + delay, event))
        self.event_list.sort(key=lambda x: x[0])
        # heappush(self.event_list,(self.t+delay,event))

    def next(self) -> Event:
        """Takes the next scheduled event from the list and advances the time accordingly.
        :return: next scheduled event
        """
        # t,e = heappop(self.event_list)
        t, e = self.event_list.pop(0)
        self.t = t
        return e

    def get_time(self) -> float:
        """Returns the current simulation time
        :return: time as a float
        """
        return self.t

    def is_empty(self) -> bool:
        """:return: true if no more events are in the list"""
        return len(self.event_list) == 0

State Variables and Parameters¶

The base class is now easily embedded to a queueing model, cosidering the three parameters $k,\mu_a,\mu_s$ and two state variables: a queue $Q$ and a number of available servers $S$ . Note how $T$ and $X$ are used for geneting simulation output.

At this point we also need to specify the distribution for $t_a$ and $t_a$ to sample new values. We leave this open by adding two provider functions as additional parameters. Both need to sample positive values with average 1. By default, the variables are sampled form an exponential distribution.

class MultipleServerQueue(States):
    def __init__(
        self,
        sim: DESSimulator,
        k: int,
        mu_a: float,
        mu_s: float,
        provider_a: Callable[[], float] | None = None,
        provider_s: Callable[[], float] | None = None,
    ):
        """Simulate a multiple server queue
        :param sim: discrete event simulator
        :param k: number of servers
        :param mu_a: average arrival time
        :param mu_s: average service time
        :param provider_a: provider for a positive arrival time (mean = 1), exp(1) by default
        :param provider_s: provider for a positive service time (mean = 1), exp(1) by default
        """
        self.sim = sim
        self.k = k
        self.mu_a = mu_a
        self.mu_s = mu_s
        self.provider_a = np.random.exponential if provider_a == None else provider_a
        self.provider_s = np.random.exponential if provider_s == None else provider_s
        self.Q = 0
        self.S = 0

    def run(self, tend: float) -> tuple[list[float], list[float]]:
        """Runs the model until a given end-time
        :param tend: end-time
        :return: list of times and queue lengths (after the event)
        """
        T = list()
        X = list()
        while not self.sim.is_empty() and self.sim.get_time() < tend:
            ev = self.sim.next()
            newevs = ev.eval(self)
            for ev, delay in newevs:  # schedule new events
                self.sim.add(ev, delay)
            T.append(self.sim.get_time())
            X.append(self.Q)
        return T, X

Events¶

In order to make the model executable, one needs to implement specific events. It quickly becomes apparent that there are exactly three events that change the behaviour of the simulation.

Arrival Event¶

On arrival of a new entity, the queue-length increases. If a server is available, the Start-Service event is instantaneously scheduled without further delay. Finally, a new Arrival event is scheduled to keep the simulation running.

Start-Service Event¶

When the service has started, an entity is taken from the queue and a server is set busy. The End-Service event is scheduled with corresponding service-time delay.

End-Service Event¶

When service has ended, a server is freed again. If entities are waiting, the next Start-Service event is scheduled.

Run Event¶

Finally, a run event needs to get the simulation started.

class ArrivalEvent(Event):
    def eval(self, sim: MultipleServerQueue):
        """Increments the queue and schedules a new arrival.
        If possible, starts service
        :param sim: object with Q and S
        :return: new events
        """
        evs = list()
        sim.Q += 1
        event = ArrivalEvent()
        delay = sim.provider_a() * sim.mu_a
        evs.append((event, delay))
        if sim.S > 0:
            event = StartServiceEvent()
            evs.append((event, 0))
        return evs


class StartServiceEvent(Event):
    def eval(self, sim: MultipleServerQueue):
        """Decrements the queue, makes the server unavailable and schedules a new end-service.
        :param sim: object with Q and S
        :return: new events
        """
        evs = list()
        sim.Q -= 1
        sim.S -= 1
        event = EndServiceEvent()
        delay = sim.provider_s() * sim.mu_s
        evs.append((event, delay))
        return evs


class EndServiceEvent(Event):
    def eval(self, sim: MultipleServerQueue):
        """Makes the server avaiable again and schedules a new start service-service if possible.
        :param sim: object with Q and S
        :return: new events
        """
        evs = list()
        sim.S += 1
        if sim.Q > 0:
            event = StartServiceEvent()
            evs.append((event, 0))
        return evs


class RunEvent(Event):
    def eval(self, sim: MultipleServerQueue):
        """Initialises the states and schedules a new arrival.
        :param sim: object with Q and S and parameter k
        :return: new events
        """
        evs = list()
        sim.Q = 0
        sim.S = sim.k
        event = ArrivalEvent()
        evs.append((event, 0))
        return evs

With the events properly defined, the first basis riund for the simulation can be executed. A plotting routine needs to be defined, to visualise the output. Executing the simulation several times also shows the impact of stochasticity on the outcome.

Note how the command np.random.seed(xxx) makes the result reproducible, regardless of randomness.

k = 2
mu_a = 1
mu_s = 2
provider_a = np.random.exponential
provider_s = np.random.exponential


def initialise_model() -> MultipleServerQueue:
    """Wrapper for initialisation of the model
    :return: instance of a multiple server queue
    """
    des = DESSimulator()  # create a DES simulator class
    msq = MultipleServerQueue(
        des, k, mu_a, mu_s, provider_a, provider_s
    )  # create MSQ instance
    msq.sim.add(RunEvent(), 0)  # add the Run event
    return msq


def plot_des_output(T, X, **kwargs) -> None:
    """Helper method to correctly display the output of the simulation (zero order hold)
    :param T: event times
    :param X: state vector after event execution
    :param kwargs: optionla kwargs passed to plt.plot
    :return:
    """
    TT = list()
    XX = list()
    for t, x in zip(T, X):
        TT.extend([t, t])
        XX.extend([x, x])
    plt.plot(TT[1:], XX[:-1], **kwargs)


# run the model once
np.random.seed(12345)
plt.figure()
msq = initialise_model()
T, X = msq.run(1000)  # run the simulation until t_end=1000
plot_des_output(T, X, color="k")  # plot
plt.title("one run")

# run the model 100 times
plt.figure()
for i in range(100):
    msq = initialise_model()
    T, X = msq.run(1000)
    plot_des_output(T, X, color="k", linewidth=0.1)
plt.title("100 runs")
plt.show()

It is evident that statistics need to be carried out in order to evaluate the results. This, however, is not trivial since all runs have a different time vector. Therefore, we need to interpolate the results

def monte_carlo_simulation(
    initialise_model: Callable[[], MultipleServerQueue], iters: int, tend: float
) -> None:
    np.random.seed(12345)
    teval = np.arange(0, tend, tend / 1000)  # evaluation times
    results = np.zeros([len(teval), iters])
    plt.figure()
    for i in range(iters):
        msq = initialise_model()
        T, X = msq.run(tend)
        xeval = np.interp(teval, T, X)  # interpolate at evaluation times
        results[:, i] = xeval
        plot_des_output(
            T,
            X,
            color="k",
            linewidth=0.1,
            zorder=0,
            label=None if i != 0 else "single run",
        )
    mean = np.mean(results, axis=1)
    std = np.std(results, axis=1)
    plt.fill_between(
        teval,
        mean - std,
        mean + std,
        color="r",
        alpha=0.3,
        zorder=1,
        label="mean +/- 1std",
    )
    plt.plot(teval, mean, color="r", zorder=1, label="mean")
    plt.legend()


k = 2
mu_a = 1
mu_s = 2
provider_a = np.random.exponential
provider_s = np.random.exponential
monte_carlo_simulation(initialise_model, iters=100, tend=1000)

Experimenting with Different Parameters¶

It stands to reason that the balance between $\mu_s/k$ and $\mu_a$ is of great importance for the stability of the system. If the number of servers $k$ is too small or the mean service time $\mu_s$ too long, the system becomes unstable and the queue grows in the long term. Otherwise, the servers are always able to clear the queue.

k = 2
mu_a = 1
mu_s = 3  # too long
provider_a = np.random.exponential
provider_s = np.random.exponential
monte_carlo_simulation(initialise_model, iters=100, tend=1000)

mu_s = 1  # short enough
monte_carlo_simulation(initialise_model, iters=100, tend=1000)

/opt/conda/lib/python3.12/site-packages/ipympl/backend_nbagg.py:392: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  self.figure.savefig(buf, format='png', dpi='figure')

Experimenting with Different Distributions¶

We furthermore explore the impact of different distributions on the outcomes, by varying the provider functions. We chose the following selection of positive distributions, each with mean equal to 1.

providers = {
    "exponential": np.random.exponential,
    "gamma": lambda: np.random.gamma(2),
    "normal": lambda: max(1 + np.random.normal() * 0.3, 0),
    "uniform": lambda: 2 * np.random.uniform(),
    "triangular": lambda: np.random.triangular(0.5, 1, 1.5),
    "deterministic": lambda: 1,
}
# plot some histograms to get an idea about how the distributions look like
plt.figure()
bins = np.arange(0, 3, 0.02)
for dist, prov in providers.items():
    vals = [prov() for i in range(10000)]
    plt.hist(vals, bins, alpha=0.3, label=dist)
plt.legend()
plt.ylim([0, 550])
plt.show()

k = 2
mu_a = 1
mu_s = 2
for dist, prov in providers.items():
    provider_a = prov
    provider_s = prov
    monte_carlo_simulation(initialise_model, iters=100, tend=1000)
    plt.title(dist)
plt.show()

In particular the skewness and variance of the distribution have a negative impact on the average queue sizes. In the optimal case of a fully deterministic simulation, the queue even remains entirely stable. This can be stated as a general property of queuing systems:

\text{if }\frac{\mu_s}{k}=\mu_a\text{, then }E(Q)<\infty\Leftrightarrow V(t_s)=V(t_a)=0.

(1)

Performance¶

Even from this simple example we easily see that the performance of the discrete event simulation is a crucial factor for the usability of the simulation - nobody wants to wait ages for termination.

The bottleneck of discret event simulators is the way how event lists are kept sorted. In above implementation, the list is re-sorted whenever a new event is added. Although, Python’s list.sort is extremely optimised and fast, it causes an effort of roughly $\mathcal{O}(N\log(N))$ whenever a new event is sorted into the list, whereas $N$ is the length of the list (likely this actually goes down to $\mathcal{O}(N)$ , since all but the new element are already in the correct order).

Anyhow, there are more efficent ways to keep things sorted. One example is the use of heap- or tree-structures, in which every element is rather a leaf on a 2-dimensional tree rather than a part of a 1-dimensional list. Python’s heapq module offers easy-to-use methods to sort new elements into a tree with effort of only $O(\log(N))$ . While this also elevates the efforts for getting the next element out of the structure from $O(1)$ (just pop from the list) to $O(\log(N))$ is provdes a considerable performance boost.

The implementation below illustrates this.

class HeapQDESSimulator(DESSimulator):
    def __init__(self):
        """Base implementation of a discrete event simulator"""
        super().__init__()
        heapify(self.event_list)

    def add(self, event: Event, delay: float):
        """Schedules a new event to the list with a given delay in the future
        :param event: event notice
        :param delay: delay greater or equal to 0
        :return:
        """
        heappush(self.event_list, (self.t + delay, event))

    def next(self) -> Event:
        """Takes the next scheduled event from the list and advances the time accordingly.
        :return: next scheduled event
        """
        t, e = heappop(self.event_list)
        self.t = t
        return e


def initialise_model_heap() -> MultipleServerQueue:
    """Wrapper for initialisation of the model with the heapq-improved des simulator
    :return: instance of a multiple server queue
    """
    des = HeapQDESSimulator()  # create a DES simulator class
    msq = MultipleServerQueue(
        des, k, mu_a, mu_s, provider_a, provider_s
    )  # create MSQ instance
    msq.sim.add(RunEvent(), 0)  # add the Run event
    return msq


k = 2
mu_a = 1
mu_s = 3  # too long
provider_a = np.random.exponential
provider_s = np.random.exponential

print("without improvement")
now = datetime.datetime.now()
for i in range(100):
    msq = initialise_model()
    msq.run(4000)
print(datetime.datetime.now() - now)

print("with improvement")
now = datetime.datetime.now()
for i in range(100):
    msq = initialise_model_heap()
    msq.run(4000)
print(datetime.datetime.now() - now)

without improvement

0:00:02.289230
with improvement

0:00:01.916892