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 holoviews as hv
import numpy as np

hv.extension("bokeh")

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 in Figure 1. The times $t_a$ and $t_s$ are hereby random variables with $E(t_a)=\mu_a$ and $E(t_s)=\mu_s$ .

Figure 1:Event graph representation of a multiple server queue.

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.

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
        :returns: 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, int, Event]] = list()  # internal event list
        heapify(self.event_list)  # turns the list into a priority queue (heap)
        self.counter = 0  # internal counter to add to the event notice to avoid
        # deadlocks

    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
        :returns:
        """
        if delay < 0:
            return
        self.counter += 1
        heappush(self.event_list, (self.t + delay, self.counter, event))

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

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

    def is_empty(self) -> bool:
        """
        :returns: 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 is None else provider_a
        self.provider_s = np.random.exponential if provider_s is 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
        :returns: 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
        :returns: 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
        :returns: 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
        :returns: 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
        :returns: 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
    :returns: 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: list[float] | list[list[float]],
    x: list | list[list],
    label: str = None,
    **kwargs,
) -> hv.Curve:
    """Helper method to correctly display the output of the simulation (zero order hold)
    :param t: event times or list of event times
    :param x: state vector or list of state vectors after event execution
    :param kwargs: optional kwargs passed to hv.Curve
    :returns: holoviews curve
    """
    if isinstance(t[0], list):  # if we only plot one DES output
        ts = list()
        xs = list()
        for ti, xi in zip(t, x, strict=True):
            ts.extend(ti)
            ts.append(np.nan)
            xs.extend(xi)
            xs.append(np.nan)
    else:
        ts = t
        xs = x
    if label is None:
        return hv.Curve((ts, xs)).opts(interpolation="steps-post", **kwargs)
    else:
        return hv.Curve((ts, xs), label=label).opts(
            interpolation="steps-post", **kwargs
        )


# run the model once
np.random.seed(12345)

msq = initialise_model()
t, x = msq.run(1000)  # run the simulation until t_end=1000
display(
    plot_des_output(t, x, color="k").opts(
        xlabel="time [s]",
        ylabel="queue length [customers]",
        xlim=(0, 1000),
        title="one run",
        responsive=True,
        height=400,
        max_width=600,
    )
)  # plot and show

# run the model 100 times
ts = list()
xs = list()
for _i in range(100):
    msq = initialise_model()
    t, x = msq.run(1000)
    ts.append(t)
    xs.append(x)
display(
    plot_des_output(ts, xs, color="k", line_width=0.5).opts(
        xlabel="time [s]",
        ylabel="queue length [customers]",
        xlim=(0, 1000),
        title="100 runs",
        responsive=True,
        height=400,
        max_width=600,
    )
)  # plot and 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
) -> hv.Overlay:
    """
    Performs and plots a Monte Carlo experiment
    :param initialise_model: method to initialise the model
    :patam iters: number of Monte Carlo samples
    :param tend: final time for the simulation
    :returns: holoviews plot
    """
    np.random.seed(12345)
    teval = np.arange(0, tend, tend / 1000)  # evaluation times
    results = np.zeros([len(teval), iters])
    ts = list()
    xs = list()
    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
        ts.append(t)
        xs.append(x)
    p1 = plot_des_output(ts, xs, color="k", line_width=0.5, label="single run")
    mean = np.mean(results, axis=1)
    std = np.std(results, axis=1)
    p2 = hv.Area(
        (teval, mean - std, mean + std),
        vdims=["lower", "upper"],
        label="mean +/- 1std",
    ).opts(color="r", alpha=0.3)

    p3 = hv.Curve((teval, mean), label="mean").opts(color="r")
    return (p1 * p2 * p3).opts(
        xlabel="time [s]",
        ylabel="queue length [customers]",
        xlim=(0, tend),
        show_legend=True,
    )


k = 2
mu_a = 1
mu_s = 2
provider_a = np.random.exponential
provider_s = np.random.exponential
display(
    monte_carlo_simulation(initialise_model, iters=100, tend=1000).opts(
        title="Monte Carlo simulation", responsive=True, height=400, max_width=600
    )
)

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
display(
    monte_carlo_simulation(initialise_model, iters=100, tend=1000).opts(
        title=f"long service time (mu_s = {mu_s})",
        responsive=True,
        height=400,
        max_width=600,
    )
)

mu_s = 1  # short enough
display(
    monte_carlo_simulation(initialise_model, iters=100, tend=1000).opts(
        title=f"short service time (mu_s = {mu_s})",
        responsive=True,
        height=400,
        max_width=600,
    )
)

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
plots = list()
bins = np.arange(0, 3, 0.02)
for _dist, prov in providers.items():
    vals = [prov() for i in range(10000)]

    hist, edges = np.histogram(vals, bins=bins, density=True)
    histogram = hv.Histogram((edges, hist), kdims="time [s]", vdims="density").opts(
        color="blue",
        alpha=0.5,
        line_color="black",
        ylim=(0, 2),
        xlim=(0, 3),
        height=300,
        responsive=True,
        max_width=500,
    )
    plots.append(histogram)
display(hv.Layout(plots).cols(3))

k = 2
mu_a = 1
mu_s = 2
for dist, prov in providers.items():
    provider_a = prov
    provider_s = prov
    display(
        monte_carlo_simulation(initialise_model, iters=100, tend=1000).opts(
            title=dist, responsive=True, height=400, max_width=600
        )
    )

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
        :returns:
        """
        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.
        :returns: 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
    :returns: 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)