Discrete-Event Simulation - Basics - Modelling and Simulation Cookbook

We give a brief introduction into the fundamendals of discrete-event simulation and give representative examples.

Objectives:

What are the components of a discrete event model and how does it work?
How is discrete event simulation fundamentally different to all other modelling strategies?
In which situations would we use this modelling concept?

from abc import abstractmethod
from heapq import heapify, heappop, heappush

import holoviews as hv
import numpy as np
import simpy
from IPython.display import HTML
from sortedcontainers import SortedList

from modelling_and_simulation_cookbook import Clock

hv.extension("bokeh")

General Idea¶

If we consider a real system or a real process which, due to its fundamental structure, can only take on discrete values – for example, a queue or a number of entities – it is clear that it cannot change its state continuously. Rather, it will always maintain its current state for a certain period and then change it abruptly. In other words, up to a certain point in time $t_{end}$ , there will have been only a finite number of times $t_1<\dots<t_k$ smaller than $t_{end}$ at which it changed its state, and, corresponding, a finite number of system states $x_0,x_1,\dots,x_k$ which the system had in between.

Most dynamic modelling approaches are based on the concept of describing by how much the states of a system change over time. Differential equation models in particular (e.g. determined using a System Dynamics approach) make this especially explicit: the derivative $dx/dt$ describes by what infinitesimal amount the system state $x$ changes as time $t$ advances by an infinitesimal amount. However, since systems with discrete states, as just described, cannot change by infinitesimally small amounts, this form of description is not necessarily the most accurate. Even temporal discretisation – i.e. updating the model in equidistant time steps as is done with e.g. difference equations – does not change improve the situation, since the points in time at which the system changes its state do not necessarily follow a predefined temporal pattern.

For this reason, discrete event simulation (DES) takes a completely different approach to describing systems: rather than asking by how much the state changes with time, DES asks when the state will next change - event-based instead of activity-based Lackner (1964). On the one hand, this form of description is sufficient because, the idea perfectly mimics how the real system updates without any approximation. In particular, the system’s time-continuous nature is preserved without any discretisation. On the other hand, this form of description is also feasible because, as already mentioned, a value-discrete system only changes at a finite number of points in time. So, knowing that there is no state-change in between, we may simply jump from one point to the next and simulate the whole system as an iterative process.

Discrete-Event Model¶

With this change of perspective, the idea behind discrete event modelling is now to describe the system in terms of those events that alter its state and the time intervals between them. As with, for example, the relationship between System Dynamics and differential equation models, a distinction must be made here between (conceptual) discrete event modelling and the discrete event simulation model (DES model or discrete event model) itself. For the prior, there are various approaches, which are determined primarily by the chosen model description language (simulation language). Notebook Discrete Event Simulation with Event Graphs gives an example for one of these languages. For now, we define what is formally meant by a DES model.

As given in Banks (2005), a DES model consists of the following components: The current state of the modelled system is descibed as a set of states. The model is updated via so-called events. Whenever an event takes place:

the event may change the state(s) of the system,
the event may schedule new events,
the event may cancel already scheduled events,
the event may change already scheduled events.

With this core definition, a conceptual DES model is already formally defined. The dynamics of the model arises naturally through the iterative processing of events, which in turn schedule new events. However, in order to better establish a link with the implementation, the definition of a DES model goes into greater depths on a number of points:

When a new event is scheduled, a so-called event notice is generated and inserted into an event list (sometimes also called event queue or future event list). The notice contains all relevant information about the event together with the scheduled time. Within the event list, the event notices are kept sorted by their scheduled time so that the first event in the list is always the one to occurs next.

In order to advance the simulation executive (sometimes called scheduler or simulation engine) takes the first event notice out of the event list, advances the simulation time (often referred to as clock) to the scheduled time from the notice (often referred to as advancing the clock), and executes the corresponding event (including state-update and creation and insorting of new event notices). Hereby, the event notice usually refers to one of a finite number of event types, which are used to cluster conceptually equivalent events - e.g. the event type increment would refer to events which increase the state variable’s value by one. After the event execution is completed, the process starts anew and is repeated until (a) the event list is empty or (b) the clock has surpassed a predefined simulation end-time. This process is visualised in Figure 1.

Figure 1:Concept of a DES model: the simulation executive takes the first event notice off the event list, updates the clock and executes the corresponding event. Hereby, states may be updated, new events may be scheduled via creation of new notices, and existing notices may be cancelled or modified.

We see that, unlike for example a differential equation model, the definition of a DES model is not given from a mathematical perspective but from a computer scientific one. In other words, the implementation of the model is already taken into account during its design. This can easily be explained historically, since the idea of discrete event modelling was always motivated by the development of computer-based simulation languages and simulation environments (see Nance (1996)). In particular, the General Purpose System Simulator (GPSS) and SIMULA were developed in the 1960s and helped to establish DES as a distinct methodology for analyzing industrial, military, and service systems. Since then, advances in computing power and software and more advanced system languages and simulators have expanded DES applications into fields including logistics, healthcare, telecommunications, and transportation.

Case-Study: Collatz Sequence, DES vs. Time-Discrete - Conceptual Model¶

In the following we give a very simple example of how a DES model is conceptualised and implemented using a case study from pure mathematics: the Collatz sequence. This sequence is defined recusively via

x_{n+1}=\left\lbrace\begin{array}{cl}x_n/2& \text{if $x_n$ is even},\\3x_n+1& \text{otherwise,}\end{array}\right.

(1)

with some initial value $x_0\in\mathbb{N}$ . It is well known from the corresponding Collatz conjecture (or $3k+1$ conjecture) which states that the sequence always eventually drops to 1 independent from its initial value. Note that, as soon as the sequence reaches 1, it will be caught in the limit cycle $1,4,2,1,\dots$ . So, this will pose our termination criterium. Note that reaching this limit cycle is subject of the still unproven famous Collatz conjecture stating, that the sequence will always eventually terminate in this cycle, independent of the initial value $x_0$ (Lothar Collatz, 1937).

For the purposes of modelling, let us assume that there is a dynamic system in which the scalar integer-valued state is updated every second according to the rules of the Collatz sequence, i.e.

\forall t\in \mathbb{R}^+, x(t):=x_{\lfloor t\rfloor}

(2)

then we can then use this as an (unusual) case-study for discrete-event modelling: The system can be modelled using a scalar integer-valued state variable $x$ and three event types:

event type	state-change	scheduling
half	$x\leftarrow \frac{x}{2}$	Depends on the value of the updated $x$ . Even and $>1$ : new event with type half in 1 second Odd and $>1$ : new event with type 3k+1 in one second.
3k+1	$x\leftarrow 3x+1$	Same as above.
run	$x\leftarrow x_0$	Same as above.

As a DES must be triggered by something, it is common practice to define initial event(s) which are manually added to the event list at the start of the simulation. In our case, the run event serves this purpose which also sets the initial value of the sequence.

At present, the conceptual model for the case study is described in a very unintuitive manner using free-form text. This situation is highly problematic in terms of producing a reproducible model description, communicating the model, and implementing it. For this reason, there are a number of simulation languages for DES (such as the event graph formalism introduced in Discrete Event Simulation with Event Graphs or the DEVS formalism) that support creating a valid and visibly appealing reproducible model desciption which can be communicated to domain experts, and allow integration into suitable DES simulators.

Simulation of DES Models¶

As we have seen, the implementation of DES models is implied by its model description. Accordingly, it is not particularly difficult to implement a DES model from scratch in any programming language – object-oriented languages such as Python or Java are the preferred choice here. However, to ensure that the self-implemented simulation executive and the event list performs well and is generic, it is advisable to give this some thought (see later). Accordingly, it is usually a better option to use a fully developed DES simulator.

Case-Study: Collatz Sequence, DES vs. Time-Discrete - Implemented Model(s)¶

Nevertheless, for the Collatz case-study we want to develop the implemented model from scratch. On the one hand, this allows us to see in detail how a simulation executive works. On the other hand, it enables us to compare the implementation with a classic time-discrete implementation of the same system (intrepreted as a difference equation), to see the structural differences in the implementation.

class CollatzModelTimeDiscrete:
    def __init__(self) -> None:
        """
        Implements the collatz sequence with iterative update
        every second (time-discrete)
        """
        pass

    def run(self, x_0: int, t_end: int = 10000) -> tuple[list, list]:
        """
        Performs the update by iteration over time
        :param x_0: starting value of the iteration
        :param t_end: maximum end-time of the simulation. If we believe in the
        conjecture, we could use inifinity here.
        :returns: two list with simulation time and state values
        """
        self.t = 0
        self.x = x_0
        ts = [self.t]
        xs = [self.x]
        for i in range(1, t_end + 1):  # loop over time
            self.t = i
            if self.x % 2 == 0:  # distinction between even and odd
                self.x = int(self.x / 2)  # cast to int just to be safe
            else:
                self.x = int(self.x * 3 + 1)
            ts.append(self.t)
            xs.append(self.x)
            if self.x == 1:
                break  # stop iteration as soon as we hit 1
        return ts, xs


class CollatzModelEventBased:
    def __init__(self):
        """
        Implements the collatz sequence with event-based update
        """
        self.event_list = list()  # create an empty event-list

    def add_event(self, event_type: str, delay: float) -> None:
        """
        Adds an event notice with the given type and specified delay (from now) to the
        event list
        :param event_type: type of the event
        :param delay: how far in the future the event should be planned
        """
        notice = (self.t + delay, event_type)
        self.event_list.append(notice)  # add the notice to the list
        # sort the list, so that the last element is the most recent one
        # note that list.pop(-1) is cheaper than list.pop(0)
        self.event_list.sort(reverse=True)

    def get_next_event(self) -> tuple[float, str]:
        """
        Returns the next event notice from the list
        :returns: event schedule time and type
        """
        return self.event_list.pop()  # get the last element from the list

    def execute_event(self, event_type: str):
        """
        Executes an event with the specified type
        :param event_type: type of the event to execute
        """
        if event_type == "half":
            self.x = int(self.x / 2)
        elif event_type == "3k+1":
            self.x = int(self.x * 3 + 1)
        elif event_type == "run":
            self.x = x_0
        else:
            raise ValueError(f"illegal event type {event_type} detected")
        # check which future event(s) should be scheduled
        if self.x > 1:
            if self.x % 2 == 0:
                self.add_event("half", 1)
            else:
                self.add_event("3k+1", 1)

    def run(self, x_0: int, t_end: int = 10000) -> tuple[list, list]:
        """
        Performs the update by jumping from update to update
        :param x_0: starting value of the iteration
        :param t_end: maximum end-time of the simulation. If we believe in the
        conjecture, we could use inifinity here.
        :returns: two list with simulation time and state values
        """
        self.event_list.clear()  # just to be safe, clear the list first
        self.t = 0
        self.x = None  # will be initialised from the run-event
        ts = list()  # can start with empty records here
        xs = list()
        self.add_event("run", 0)  # add the initial event to the list
        while self.t < t_end and len(self.event_list) > 0:
            time, event_type = self.get_next_event()
            self.t = time
            self.execute_event(event_type)
            ts.append(self.t)
            xs.append(self.x)
        return ts, xs


mdl_1 = CollatzModelTimeDiscrete()
mdl_2 = CollatzModelEventBased()
x0 = 27
for x_0 in [10, 27, 80]:
    [tt, xx] = mdl_1.run(x_0)
    c1 = hv.Curve((tt, xx)).opts(
        xlabel="time",
        ylabel="value of x",
        responsive=True,
        max_width=400,
        height=400,
        interpolation="steps-post",
        title=f"time-discrete model {x_0=}",
    )
    [tt, xx] = mdl_2.run(x_0)
    c2 = hv.Curve((tt, xx)).opts(
        xlabel="time",
        ylabel="value of x",
        responsive=True,
        max_width=400,
        height=400,
        interpolation="steps-post",
        title=f"DES model {x_0=}",
    )
    display(c1 + c2)

The graphs of the results show immediatly why this mathematical problem is so extraordinarily fascinating. The number of steps it takes for the sequence to drop to 1 seems to depend almost arbitrarily of the chosen initial value and is virtually impossible to predict. Although the number increases on average with to the size of the initial value, the volatility is astonishing. Note how the sequence starting in 27 surpasses values of 9000 before eventually dropping to the 1-4-2 cycle.

However, the key takeaway from the case study is, of course, the implementation as a DES model. One immediately recognises that a key distinguishing feature of the DES model is that the update loop does not run over the simulation time, but rather iterates through the event list. It is quite common to use some while loop here, as the number of events up to the simulation end-time is usually unpredictable.

Moreover, note how the interpolation method interpolation='steps-post' was used to make sure that the plots are displayed with a zero-hold and not by linear spline interpolation between the states. This is important for discrete event simulation, since the states are defined to remain constant between any two consecutive events.

Collatz Sequence, DES vs. Time-Discrete - Flexibility of DES¶

While the DES strategy clearly produces overheads in this particular case compared to the time-discrete strategy, it gives a lot of freedom for designing the update - not only in terms of how the update evaluates but also when it happens. For example, one quickly notices that the update $3x+1$ for an odd $x$ will always result in an even number and lead to consequent halving of the number. In other words:

x_{n+2}=\frac{3x_n+1}{2}

(3)

if $x_n$ is odd. In case we only want to find the number of iterations until the sequence collapses, we may skip computing $x_{n+1}$ in these cases. This feature is comparably tricky to include in a sole time-discrete approach, because it would be necessary to skip selected time-steps. However, as seen below, it is a trivial change in the event-based architecture.

class CollatzModelEventBased2(CollatzModelEventBased):
    def execute_event(self, event_type: str):
        """
        Executes an event with the specified type
        :param event_type: type of the event to execute
        """
        if event_type == "half":
            self.x = int(self.x / 2)
        elif event_type == "(3k+1)/2":
            self.x = int((self.x * 3 + 1) / 2)  # changing from 3x+1 to (3x+1)/2
        elif event_type == "run":
            self.x = x_0
        else:
            raise ValueError(f"illegal event type {event_type} detected")
        # check which future event(s) should be scheduled
        if self.x > 1:
            if self.x % 2 == 0:
                self.add_event("half", 1)
            else:
                self.add_event(
                    "(3k+1)/2", 2
                )  # changing from a delay of 1 to a delay of 2


mdl_1 = CollatzModelEventBased()
mdl_2 = CollatzModelEventBased2()
x_0 = 27
[tt, xx] = mdl_1.run(x_0)
c1 = hv.Curve((tt, xx)).opts(
    xlabel="time",
    ylabel="value of x",
    responsive=True,
    max_width=400,
    height=400,
    interpolation="steps-post",
    title=f"DES model {x_0=}",
)
[tt, xx] = mdl_2.run(x_0)
c2 = hv.Curve((tt, xx)).opts(
    xlabel="time",
    ylabel="value of x",
    responsive=True,
    max_width=400,
    height=400,
    interpolation="steps-post",
    title=f"(3k+1)/2 DES model {x_0=}",
)
display(c1 + c2)
nums1 = list()
nums2 = list()
for x_0 in range(1, 1000):
    nums1.append(mdl_1.run(x_0)[0][-1])
    nums2.append(mdl_2.run(x_0)[0][-1])
c1 = hv.Scatter([(i + 1, x) for i, x in enumerate(nums1)]).opts(
    xlabel="x_0",
    ylabel="steps to reach 1",
    responsive=True,
    max_width=400,
    height=400,
    title="DES model",
)
c2 = hv.Scatter([(i + 1, x) for i, x in enumerate(nums2)]).opts(
    xlabel="x_0",
    ylabel="steps to reach 1",
    responsive=True,
    max_width=400,
    height=400,
    title="(3k+1)/2 DES model",
)
display(c1 + c2)

While the model shows different dynamics due to skipping the $x_{n+1}=3x_n+1$ values, the number of steps for dropping to 1 are identical. The lower plots show the fascinated dynamics between the starting value $x_0$ and the time until the model state reaches 1 for the first time.

Although the Collatz sequence is not a typical example of when DES is used for modelling or simulation, even in this simple case-study one can already see the advantages of not only defining how an update takes place, but also when it should occur.

Stochasticity in DES¶

Stochasticity is one of the most important components in DES models. Although, as the Collatz example shows, it is possible to define DES models deterministically, the vast majority of DES models are stochastic. It should be noted here that stochasticity may relate not only to the updating of the state, but also, and above all, to the time intervals between updates, i.e. the scheduling times of the events.

Collatz Sequence, DES vs. Time-Discrete - Stochasticity¶

Altough the original Collatz problem is deterministic, there are indeed stochastic variants of it which we could embed in our case-study, for example, the one published by Ingo Althofer on his homepage where the update for odd numbers is either $3x_n+1$ or $3x_n-1$ with equal probability. However, this randomness relates purely to the update of the state and could be depicted with a time discrete model as well. In our version, however, we want to modify the times of the state-updates to incorporate stochasticity to better see the advantages of the approach.

Suppose, the investigated system $x$ does not update reliably every second, but updates my occur at any time with a certain constant likelihood. A system like this would then be modelled using a so-called Markov process with a certain rate - the higher the rate, the higher the probability that an update occurs within a given interval. Theory about these processes tells us, that the time between any two updates of a Markov proess is distributed exponentially $Exp(\lambda)$ . Hereby the parameter $\lambda$ of the exponential disribution is equal to the transition rate of the Markov process and the expected number of updates per time.

We use this information to modify our model. Instead of scheduling an update event every second instead schedule it with a time delay of $X\sim E(1)$ seconds.

class CollatzModelEventBasedStochastic(CollatzModelEventBased):
    def execute_event(self, event_type: str):
        """
        Executes an event with the specified type
        :param event_type: type of the event to execute
        """
        if event_type == "half":
            self.x = int(self.x / 2)
        elif event_type == "3k+1":
            self.x = int(self.x * 3 + 1)
        elif event_type == "run":
            self.x = x_0
        else:
            raise ValueError(f"illegal event type {event_type} detected")
        # check which future event(s) should be scheduled
        if self.x > 1:
            if self.x % 2 == 0:
                self.add_event(
                    "half", np.random.exponential(1)
                )  # changing from a delay from 1 to Exp(1)
            else:
                self.add_event(
                    "3k+1", np.random.exponential(1)
                )  # changing from a delay from 1 to Exp(1)


np.random.seed(12345)
mdl_1 = CollatzModelEventBased()
mdl_2 = CollatzModelEventBasedStochastic()
x_0 = 12
[tt, xx] = mdl_1.run(x_0)
c1 = hv.Curve((tt, xx)).opts(
    xlabel="time",
    ylabel="value of x",
    responsive=True,
    max_width=400,
    height=400,
    interpolation="steps-post",
    title=f"DES model {x_0=}",
)
[tt, xx] = mdl_2.run(x_0)
[tt2, xx2] = mdl_2.run(x_0)
c2 = (
    hv.Curve((tt, xx)).opts(interpolation="steps-post")
    * hv.Curve((tt2, xx2)).opts(interpolation="steps-post")
).opts(
    xlabel="time",
    ylabel="value of x",
    responsive=True,
    max_width=400,
    height=400,
    title=f"Stochastic DES model (2 runs) {x_0=}",
)
display((c1 + c2).opts(shared_axes=False))

While the overall state updates are identical, the stochasticity of the schedule times causes temporarily streched and shifted curves. Moreover, every execution of the stochastic model leads different results due to the use of randomness. We also see the impact of the skewness of the exponential distribution, i.e. that short timespans are quite likely, but very long timespans are also possible.

Computational Effort of DES¶

Above case study also raises interesting questions related to the comparison of computational efforts.

In a time-discrete update, the computational costs of the approach always increase linearly with the simulation time $t_{end}$ divided by the length $h$ of the used time-steps. In the Collatz example above, $h=1$ , however thinking about numerical simulation of differential equation models (e.g. using the Euler method) choosing a small stepsize $h$ is important for getting a proper approximation to the analytic solution. Summarising,

costs=\mathcal{O}\left(\frac{t_{end}}{h}\right)

(4)

for implementation of time discrete(tised) models.

In DES approaches, we do not use any time-steps but iterate of the event-list. Let $\mathbb{E}(h)$ refer to the average time-distance between any two consecutive events, then we get, again, a computational effort of $\mathcal{O}(t_{end}/\mathbb{E}(h))$ . However, in this computation we are neglecting a potentially critical contributer to computation time: efforts for magagement of the event list. On every update step, the simulation executive must (a) take and remove the next-scheduled event notice from the event list, (b) insort any new created event notices, and (c) search and modify/delete already scheduled event notices. The computational effort of theses processes varies with the used data structure, however, it can be generally stated that at least one of them has a computational effort which depends on the average length $\mathbb{E}(K)$ of the event list. So, generally speaking,

costs=\mathcal{O}\left(\frac{t_{end}\cdot f(\mathbb{E}(K))}{\mathbb{E}(h)}\right),

(5)

for some strict monotonically increasing function $f$ . The Collatz example was extreme in this sense, since $K\equiv 1$ all the time. However, in the more general case, event lists can become really long in models with lots of different event types, in particular, in agent-based models with event-based update. On the one hand, one should consider this when opting for or against an event-based update strategy compared to a discret(ised) one in case both options are viable. On the other hand, one should have a proper look at how event lists and the corresponding adding and removing of event notices is implemented in the corresponding software (i.e. what is “ $f$ ”).

We give some examples for the efforts of handling event lists with $K$ elements with different data structures / methods in Python:

data structure for event list	example Python object (alternatives)	find and remove next event notice	insert new event notice at the right place	overall $f$
list (sorted) or linked list	`list` (`deque`, `np.ndarray`)	`pop` $\mathcal{O}(1)$	`append` $\mathcal{O}(1)$ + `sort` $\mathcal{O}(K)$	$\mathcal{O}(K)$
list (unsorted)	`list` (`deque`, `np.ndarray`)	`min` $\mathcal{O}(K)$ + `remove` $\mathcal{O}(K)$	`append` $\mathcal{O}(1)$	$\mathcal{O}(K)$
set	`set`	`find` $\mathcal{O}(K)$ + `remove` $\mathcal{O}(1)$	`add` $\mathcal{O}(1)$	$\mathcal{O}(K)$
heap	`heapq` (`queue.PriorityQueue`)	`heappop` $\mathcal{O}(ln(K))$	`heappush` $\mathcal{O}(ln(K))$	$\mathcal{O}(ln(K))$
binary tree	`sortedcontainers.SortedList`	`pop` $\mathcal{O}(ln(K))$	`add` $\mathcal{O}(ln(K))$	$\mathcal{O}(ln(K))$

However, in a Jupyter Notebook, we don’t necessarily have to rely solely on theoretical considerations and embed some actual benchmarks...

class DESSimulator:
    """
    Abstract class for applying a standartized test
    """

    @abstractmethod
    def reset(self) -> None:
        pass

    @abstractmethod
    def add(self, value: float) -> None:
        pass

    @abstractmethod
    def get(self) -> float:
        return 0.0


class List1Simulator(DESSimulator):
    """
    Uses list-sort and pop
    """

    def __init__(self):
        self.event_list = list()

    def reset(self):
        self.event_list.clear()

    def add(self, value: float):
        self.event_list.append(value)
        self.event_list.sort(reverse=True)

    def get(self):
        return self.event_list.pop()


class List2Simulator(DESSimulator):
    """
    Uses an unsorted list, searches for the most recent element and removes it
    """

    def __init__(self):
        self.event_list = list()

    def reset(self):
        self.event_list.clear()

    def add(self, value: float):
        self.event_list.append(value)

    def get(self):
        min_value = min(self.event_list)
        self.event_list.remove(min_value)
        return min_value


class SetSimulator(DESSimulator):
    """
    Same as List2Simulator, but with a set object for cheaper removal
    """

    def __init__(self):
        self.event_list = set()

    def reset(self):
        self.event_list.clear()

    def add(self, value: float):
        self.event_list.add(value)

    def get(self):
        min_value = min(self.event_list)
        self.event_list.remove(min_value)
        return min_value


class HeapSimulator(DESSimulator):
    """
    Uses a heap structure (heapify) and push and pop
    """

    def __init__(self):
        self.event_list = list()
        heapify(self.event_list)

    def reset(self):
        self.event_list = list()
        heapify(self.event_list)

    def add(self, value: float):
        heappush(self.event_list, value)

    def get(self):
        return heappop(self.event_list)


class TreeSimulator(DESSimulator):
    """
    Uses a tree structure (SortedList) with add and pop
    """

    def __init__(self):
        self.event_list = SortedList(
            key=lambda x: -x
        )  # add key to make sure that the list is sorted in reverse

    def reset(self):
        self.event_list.clear()

    def add(self, value: float):
        self.event_list.add(value)

    def get(self):
        return self.event_list.pop()


def test(
    des_simulator: DESSimulator, queue_size: int, iters: int, iters_2: int = 200
) -> tuple[float, float]:
    """
    Evaluates a time-test for the simulator.
    1) the event list is filled with `queue_size` events
    2) a total of `iters_2` event notices are added to the event list
    3) the next `iters_2` event notices are pulled from the list
    4) steps 2 and 3 are repeated `iters` iterations
    Times for adding and pulling are measured and returned separately
    :param des_simulator: instance of the simulator to test.
        Must implement `DESSimulator`
    :param queue_size: initial number of events to put into the event list
        (the higher, the more problematic is O(N) and O(ln(N)) effort)
    :param iters: number of repetitions of the filling/emptying process
        (the higher, the more accurate is the time measurement)
    :param iters_2: number of events added/removed from the queue every iteration
        (the lower, the more inaccurate the time measurement,
        should be well below queue_size though)
    :returns: time (in seconds) used for adding and removing of events, respectively
    """
    np.random.seed(12345)
    des_simulator.reset()
    for _i in range(queue_size):
        des_simulator.add(np.random.random())
    c = Clock()
    t_add = 0.0
    t_get = 0.0
    for _i in range(iters):
        c.start()
        for _j in range(iters_2):
            des_simulator.add(np.random.random())
        t_add += c.stop(show_message=False).total_seconds() * 1000
        c.start()
        for _j in range(iters_2):
            des_simulator.get()
        t_get += c.stop(show_message=False).total_seconds() * 1000
    return t_add, t_get


def test_order(
    des_simulator: DESSimulator, queue_sizes: list[int], iters: int, iters_2: int = 200
) -> str:
    """
    iteratively calls `test` with different queue sizes to determine how the
    runtime of the operations depend on N, i.e. O(1),O(ln(N)) or O(N)
    :param des_simulator: instance of the simulator to test.
        Must implement `DESSimulator`
    :param queue_sizes: list of event list sizes to test.
        Should span different magnitudes to properly evaluate the order
    :param iters: number of repetitions of the filling/emptying process
        (the higher, the more accurate is the time measurement)
    :param iters_2: number of events added/removed from the queue every iteration
        (the lower, the more inaccurate the time measurement,
        should be well below queue_size though)
    :returns: html-formatted string with the table of the test results
    """
    t_adds = np.zeros([len(queue_sizes), 3])
    t_gets = np.zeros([len(queue_sizes), 3])
    norm_queue = 1000
    for i, q in enumerate(queue_sizes):
        t_add, t_get = test(des_simulator, q, iters, iters_2)
        t_adds[i, :] = [
            t_add,
            t_add * np.log(norm_queue) / np.log(q),
            t_add * norm_queue / q,
        ]
        t_gets[i, :] = [
            t_get,
            t_get * np.log(norm_queue) / np.log(q),
            t_get * norm_queue / q,
        ]
    # assumption: if the process scales with order f(N), then the runtime t_N
    # divided by f(N) should be constant
    # so we compute the coefficient of variation and assume that the "f" with the
    # smallest coefficient is the most likely complexity order
    var_coeff = np.std(t_adds, axis=0) / np.mean(t_adds, axis=0)
    add_mins = np.where(np.logical_or(var_coeff == min(var_coeff), var_coeff < 0.2))
    var_coeff = np.std(t_gets, axis=0) / np.mean(t_gets, axis=0)
    get_mins = np.where(np.logical_or(var_coeff == min(var_coeff), var_coeff < 0.2))
    html_table = f"<b>{type(des_simulator).__name__}</b>"
    html_table += "<table><tr><td>queue-length</td><td><b>add time</b> [ms]</td>"
    html_table += f"<td>[ms*log({norm_queue})/log(K)]</td>"
    html_table += f"<td>[ms*{norm_queue}/K]</td><td><b>get time</b> [ms]</td>"
    html_table += f"<td>[ms*log({norm_queue})/log(K)]</td>"
    html_table += f"<td>[ms*{norm_queue}/K]</td></tr>"
    for i, q in enumerate(queue_sizes):
        html_table += f"<tr><td>{q}</td>"
        for j in range(3):
            col = "red" if j in add_mins[0] else "black"
            html_table += f'<td style="color:{col}">{t_adds[i][j]:0.1f}</td>'
        for j in range(3):
            col = "red" if j in get_mins[0] else "black"
            html_table += f'<td style="color:{col}">{t_gets[i][j]:0.1f}</td>'
        html_table += "</tr>"
    html_table += "</table>"
    return html_table


iters = 100
queue_sizes = [1000, 5000, 10000]
display(HTML(test_order(List1Simulator(), queue_sizes, iters)))
display(HTML(test_order(List2Simulator(), queue_sizes, iters)))
display(HTML(test_order(SetSimulator(), queue_sizes, iters)))
queue_sizes = [1000, 10000, 50000, 100000]
display(HTML(test_order(HeapSimulator(), queue_sizes, iters)))
display(HTML(test_order(TreeSimulator(), queue_sizes, iters)))

The test displays times for adding and getting events from corresponding event lists, absolute and relative to the logarithm and the total size of the queue. The red column(s) in the tables have, according to the benchmark results, the highest likelihood to represent the order of the computational effort. The methods with tree or heap structure are clearly superior to the other approaches in terms of computational performance. In both cases, the computational effort increases, at most, with the logarithm of the queue size. For the list and set approaches, at least one of the two operations scales linearly with the queue length.

When it comes to searching and modification/removal of events, basically all options are equally good/bad, since $\mathcal{O}(N)$ cannot be avoided. As a result, it is always a good idea to avoid these in the model design, when possible.

DES Simulators¶

As previously said, one should carefully consider how to implement a DES model in order to avoid performance bottlenecks. Note that due to the common stochasticity of DES models,a single simulation run is often of very limited informative value - so high computational performance is even more important to create mulitple runs (Monte Carlo simulation). Anyway, using an existing simulator is usually a good strategy.

In most cases, simulators are specifically tailored to the chosen simulation / model description language. For example, the SIGMA simulator supports the EventGraph formalism, PowerDevs supports the DEVS formalism, and AnyLogic supports a formalism based on process diagrams.

In Python, the simulator SimPy is likely the most frequently used. It is a general purpose simulation environment, which offers a multitude of methods for discrete event simulation without bein explicitly specificied for amy particular model description language.

class CollatzModelSimpyStochastic:
    def __init__(self) -> None:
        """
        Performs the collatz sequence with event-based update simulated with SimPy
        """
        pass  # no need to do anything here

    def collatz_process(self, x_0: int):
        """
        Actual process implementing the Collatz sequence.
        :param x_0: initial value of the sequence
        """
        self.x = x_0

        self.ts.append(self.env.now)
        self.xs.append(self.x)

        while self.x > 1:
            yield self.env.timeout(
                np.random.exponential(1)
            )  # this creates the time-delay between any two events
            if self.x % 2 == 0:
                self.x = int(self.x / 2)
            else:
                self.x = int(3 * self.x + 1)
            self.ts.append(self.env.now)
            self.xs.append(self.x)

    def run(self, x_0: int, t_end=10000):
        """
        Performs the update by jumping from update to update
        :param x_0: starting value of the iteration
        :param t_end: maximum end-time of the simulation. If we believe in the
        conjecture, we could use inifinity here.
        :returns: two list with simulation time and state values
        """
        self.env = (
            simpy.Environment()
        )  # create the SimPy environment including event scheduled
        self.x = None  # setup the state
        self.ts = list()  # setup the output
        self.xs = list()  # setup the output
        self.env.process(self.collatz_process(x_0))  # assign the model process
        self.env.run(until=t_end)  # perform the simulation
        return self.ts, self.xs


np.random.seed(12345)
mdl_1 = CollatzModelEventBasedStochastic()
x_0 = 12
[tt, xx] = mdl_1.run(x_0)
[tt2, xx2] = mdl_1.run(x_0)
c1 = (
    hv.Curve((tt, xx)).opts(interpolation="steps-post")
    * hv.Curve((tt2, xx2)).opts(interpolation="steps-post")
).opts(
    xlabel="time",
    ylabel="value of x",
    responsive=True,
    max_width=400,
    height=400,
    title=f"Stochastic DES model (2 runs) {x_0=}",
)
np.random.seed(12345)
mdl_2 = CollatzModelSimpyStochastic()
[tt, xx] = mdl_2.run(x_0)
[tt2, xx2] = mdl_2.run(x_0)
c2 = (
    hv.Curve((tt, xx)).opts(interpolation="steps-post")
    * hv.Curve((tt2, xx2)).opts(interpolation="steps-post")
).opts(
    xlabel="time",
    ylabel="value of x",
    responsive=True,
    max_width=400,
    height=400,
    title=f"Stochastic DES model with SimPy (2 runs) {x_0=}",
)
display((c1 + c2).opts(shared_axes=False))

Usage of DES¶

Given the features discussed, DES is particularly advantageous for systems that change only at specific points in time, however it should be noted here that DES can only truly shine with its full potential if these points do not follow a predefined time cycle and/or include stochasticity. Otherwise, a time-discrete approach in which a time-loop steers the simulation is usually more suitable.

Queueing systems are undoubtedly the prime example of the application of DES. The state (the length of the queue) is discrete and changes only when an entity arrives or is served (see Figure 2). We refer to Discrete Event Simulation with Event Graphs for more queueing case examples.

Figure 2:Parade example of a system which is best modelled using DES: a queue with server.

While systems with discrete values are best suitet for DES it should be noted that also value-continuous systems can be modelled with this strategy, when interpreted in the right way. A good example for this is a bouncing ball (see Figure 3, panel [a]): the moments at which the ball touches the ground could be modelled as events, at which the ball’s velocity changes its sign and is damped. The time interval until the next impact (event) can then be derived quite easily from the laws of physics. In other words, DES is applicable when a continuous variable changes only in discrete steps (or when the discrete change is the essential element of the modelling). Another example where DES is applicable for value-continuous systems is given by a train (Figure 3, panel [b]) which, although it is generally always located at a continuous position, is reduced to the states: ‘at station A’, ‘between A and B’ and ‘at station B’. In other words, DES is applicable for value-continuous systems, when it is reasonable to discretise its value.

Figure 3:Examples for systems with continuous values which are still eligble for DES modelling: a bouncing ball modelled via bounce-events (lefT) and a train’s position simplified to “in station A/B” and “in between stations” modelled with departure and arrival events.

References¶

Lackner, M. R. (1964). Digital Simulation and System Theory [Techreport]. System Development Corporation, SDC SP-1612.
Banks, J. (2005). Discrete event system simulation. Pearson Education India.
Nance, R. E. (1996). A history of discrete event simulation programming languages. In History of Programming Languages—II (pp. 369–427). Association for Computing Machinery. 10.1145/234286.1057822