Simulation Circle - Modelling and Simulation Cookbook

In this chapter we present the simulation circle and give an explorative example.

Objectives:

How are the steps involved in a classic modelling and simulation study?
Why is this considered a “loop”?
What is validation and verification and what is the difference between the two?

import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import solve_ivp
from scipy.interpolate import interp1d

Simulation Circle Concept¶

The simulation circle, sometimes also referred to as the modelling circle, modelling/simulation lifecycle or modelling/simulation cycle, depicts the generalised history of a model’s development from the research question to the simulation experiments in a diagram. Although this usually includes the steps already described in the Basics chapter, the most important feature of this diagram is the concept of the cycle, which is intended to illustrate that the development of a model is not a single pipeline, but an iterative process. There are many versions of this modelling cycle, with different versions having different focuses. For historical and didactical reasons, this course will follow the diagram below. It dates back, at least, to the 1990s and has eversince been used for modelling and simulation education at TU Wien. Its original source is unclear but heavily tied to Felix Breitenecker, who taught corresonding courses in the past. Finally, the diagram has been edited by Martin Bicher in the 2020s. In its core it matches other noteworthy examples like those in Banks (2005), Figure 1.3, or Sargent (2010), Figure 2. Of the two, the prior lays additional emphasis on the integration of data, whereas the latter puts a particular focus on the strict differentiation between conceptual and implemented model.

In the following we demonstrate the steps of the simulation circle based on a simple example from population dynamics.

Problem¶

Every model must follow a certain pragmatism, i.e. a certain purpose. In modelling and simulation studies, this modelling purpose is given by a problem, the study should help to solve. In most cases, this problem is exogenous and comes from an expert in the respective application domain. As a result, usually the modeller must first analyse the problem and break it down into specific, measurable questions that can be solved using modelling and simulation. We give a simple fictional example from agriculture.

Case Study: Predator-Prey - Problem¶

A farmer hires us as experienced modellers to help him with a problem with his harvest. An invasive rodent species is increasingly causing a decline in harvests as the animals feed on the roots. The farmer sees the use of chemicals that reduce the fertility of the rodents as last chance. The farmer’s question is now: does the use of chemicals make sense and when should they be used?

Modelling¶

The next step is to design the model. As mentioned above, the modeller decides on a modelling approach and follows it through to the conceptual model. The decision on the approach depends on many factors, but above all on the structure of the system and the given research problem. However, aspects as available data, computational resources, specific preferences and flexibility requirements must not be underestimated.

Case Study: Predator-Prey - Modelling¶

In addressing this question, the modeller follows a System Dynamics approach. We will not go into detail here about how the approach works in detail, but refer to the corresponding chapter instead. The involved processes are roughly the following:

The modeller learns from domain experts (farmers, foresters, hunters, etc.) that the population dynamics of rodents are primarily regulated by the regional fox population. A dataset (see below) obtained from counts of local hunters furthermore suggest an ocillating behaviour of the corresponding population size.
He therefore chooses a macroscopic modelling approach with two main system components: the rodent and the fox population.
Through system analysis, the modeller discovers various system components that are considered relevant to the model. The modeller depicts these in the form of the Causal Loop Diagram depicted below (left).
The modeller analyses the role of the components in the system and identifies stocks, flows, links, auxiliaries and parameters, and puts them together as the Stock and Flow Diagram depicted below (right).
The Stock and Flow diagram can be translated to a differential equation model with two coupled differential equations. Let $x(t)$ and $y(t)$ stand for the number of prey and predator animals (i.e. the values of the two stocks), then the system can be described as
$\frac{dx}{dt}(t)=x(t)(\alpha_1\cdot(1-\xi(t))-\beta_1)-x(t)y(t)\gamma,$
(1)

$\frac{dy}{dt}(t)=y(t)(\alpha_2-\beta_2)+x(t)y(t)\gamma\delta.$
(2)
Hereby, $\alpha_1,\beta_1$ refers to the birth- and death-rate of the prey ([births/(prey*month)] and [deaths/(prey*month)]), $\alpha_2,\beta_2$ to the corresponding rates of the predators. Moreover, $\gamma$ stands for the hunting-rate [deaths/(prey*predator*month)], and $\delta$ stands for the impact of a hunted prey animal on the predators, i.e. how many new predators are born per hunted prey animal [births/death]. Finally, $\xi(t)$ stands for the multiplicative impact of the chemicals on the birth rate. It is a function of time which is 0 outside of the intervention timeframe, and a positive value $\xi$ during reducing the rate by $100\xi\%$ . The additional specification $x(t)=x_0$ and $y(0)=y_0$ finally results in an initial-value problem with unique solution.

Note, that without considering the intervention, the system matches the famous predator-prey differential equations published by Lotka and Volterra in 1926 which are also called “Lotka-Volterra equations”.

As mentioned before, the following dataset is available, providing estimates of the size of the fox population over the last 60 months. Our goal is to use this data for model parametrisation.

t_data = [
    0.0,
    4.0,
    8.0,
    12.0,
    16.0,
    20.0,
    24.0,
    28.0,
    32.0,
    36.0,
    40.0,
    44.0,
    48.0,
    52.0,
    56.0,
    60.0,
]
pred_data = [
    x * 1000
    for x in [
        1.0,
        2.4,
        0.8,
        0.2,
        0.1,
        0.2,
        2.0,
        1.6,
        0.4,
        0.2,
        0.1,
        1.8,
        1.6,
        0.4,
        0.2,
        0.1,
    ]
]


def create_figure():
    """Just creates a new figure with suitable size"""
    plt.figure(figsize=(10, 5))


def plot_pred_data(**kwargs):
    """Plots the predator datapoints.
    :param kwargs: Standard keywords for pypot.plot.
    If not specified, this method will plot black points
    :return:
    """
    if "label" not in kwargs:
        kwargs["label"] = "foxes (counts)"
    if "linestyle" not in kwargs:
        kwargs["linestyle"] = ""
    if "marker" not in kwargs:
        kwargs["marker"] = "o"
    if "color" not in kwargs:
        kwargs["color"] = [0, 0, 0]
    plt.plot(t_data, pred_data, **kwargs)


def finish_plot():
    """Convenience method to create a legend and name the axis"""
    plt.legend()
    plt.xlabel("months")
    plt.ylabel("individuals")


create_figure()
plot_pred_data()
finish_plot()

Model Analysis¶

Before implementing the model in a programming environment, it is important to first perform a formal analysis of the conceptual model. This not only helps to validate the model, but also allows generalisable statements about the system to be made that are not based on simulations.

Case Study: Predator-Prey - Model Analysis¶

Looking at the differential equations already allows us to make some very fundamental observations. With $y_0=0$ also $y(t)\equiv 0$ , same for $x_0=0$ and $x(t)\equiv 0$ , since the respective right-hand sides disappear. That means, if a population was absent at the beginning, it remains absent. In direct consequence of the disappearance of the interaction terms $x(t)y(t)$ , the other variable is left on its own and is (neglecting the intervention) determined by

\frac{dx}{dt}(t)=x(t)(\alpha_1-\beta_1),

(3)

and

\frac{dy}{dt}(t)=y(t)(\alpha_2-\beta_2),

(4)

respectively. This linear differential equation will either lead to exponential growth, when the birth-rate exceed the death-rate, or to decay, in the other case. Since we may expect that prey animals thrive in the absence of a predator, we may assume $\alpha_1>\beta_1$ . The predators, however, are expected to die out if the corresponing prey is not present, i.e. we may assume $\alpha_2<\beta_2$ .

Numerics / Programming¶

In the next step, the modeller chooses the programming environment in which the model is to be implemented. A basic implementation is then carried out. This involves making additional decisions. On the one hand, initial parameter values must be defined. On the other hand, it may be necessary to select and properly configure numerical approximation algorithms.

Case-Study: Predator-Prey (Numerics / Programming)¶

In the specific application, the modeller can either implement the Stock and Flow Diagram directly in a System Dynamics simulator or use the differential equation representation for direct implementation in a programming language which supports their numerical solutions. An example for the prior can be seen in the image below, which shows a screenshot of the model implemented in the AnyLogic simulator.

The class below, however, represents a direct implementation of the differential equation system using Python. For numerical solution, the class uses a Runge-Kutta method (4th order with stepsize control by Dormand and Prince).

class PredPreySim:
    def __init__(
        self,
        alpha_1: float,
        beta_1: float,
        alpha_2: float,
        beta_2: float,
        gamma: float,
        delta: float,
        xi: float = 0.0,
        xi_time: [float, float] = None,
    ):
        """Initialises a new predator-prey simulation model
        :param alpha_1: birth rate prey
        :param beta_1: death rate prey
        :param alpha_2: birth rate predators
        :param beta_2: death rate predators
        :param gamma: number of prey anmials eaten per prey and per predator animal
        :param delta: number of newborn predator from one successful prey hunt
        :param xi: effectiveness of the chemical intervention on the prey's birth rate
        :param xi_time: timeframe in which the intervention is active
        """
        self.alpha_1 = alpha_1
        self.beta_1 = beta_1
        self.alpha_2 = alpha_2
        self.beta_2 = beta_2
        self.gamma = gamma
        self.delta = delta
        self.xi = xi
        if xi_time is None:
            self.xi_time = [-1, -1]
        else:
            self.xi_time = xi_time

    def rhs(self, t: float, xy: np.ndarray) -> np.ndarray:
        """Right hand side of the differential equation
        :param t: current time
        :param xy: current number of prey and predator animals as vector
        :return: vector of change of the prey and predator animals
        """
        prey = xy[0]
        pred = xy[1]
        eaten = pred * prey * self.gamma
        newpred = eaten * self.delta
        if self.xi_time[0] < t <= self.xi_time[1]:
            preyBirth = prey * self.alpha_1 * (1 - self.xi)
        else:
            preyBirth = prey * self.alpha_1
        preyDeath = prey * self.beta_1 + eaten

        predBirth = pred * self.alpha_2 + newpred
        predDeath = pred * self.beta_2
        return np.array([preyBirth - preyDeath, predBirth - predDeath])

    def run(
        self, tend: float, prey_0: float, pred_0: float
    ) -> [np.ndarray, np.ndarray, np.ndarray]:
        """Runs the simulation for a given number of time units and initial population figures
        :param tend: end time of the simulation
        :param prey_0: initial prey count
        :param pred_0: initial predtor count
        :return: three vectors with equal length: reference time points and time series' of the prey and predator animals
        """
        t_model = np.arange(0, tend + 0.1, 0.1)
        sol_model = solve_ivp(
            self.rhs, [0, tend], np.array([prey_0, pred_0]), t_eval=t_model
        )
        return sol_model.t, sol_model.y[0, :], sol_model.y[1, :]


def plot_pred_model(t_model: np.ndarray, pred_model: np.ndarray, **kwargs):
    """Convenience method to plot the time-series of the predator population from a predator-prey model
    :param t_model: reference time points
    :param pred_model: time series of the predator animals
    :param kwargs: stadard kwarguments of pyplot.plot. By default, uses red solid lines
    """
    if "label" not in kwargs:
        kwargs["label"] = "foxes (simulation)"
    if "color" not in kwargs:
        kwargs["color"] = [1, 0, 0]
    plt.plot(t_model, pred_model, **kwargs)


def plot_prey_model(t_model: np.ndarray, prey_model: np.ndarray, **kwargs):
    """Convenience method to plot the time-series of the prey population from a predator-prey model
    :param t_model: reference time points
    :param pred_model: time series of the prey animals
    :param kwargs: stadard kwarguments of pyplot.plot. By default, uses blue solid lines
    """
    if "label" not in kwargs:
        kwargs["label"] = "rodents (simulation)"
    if "color" not in kwargs:
        kwargs["color"] = [0, 0, 1]
    plt.plot(t_model, prey_model, **kwargs)

Basic Simulation and Simulation Results¶

With reasonably estimated parameters, basic simulations are carried out in the next step. These are intended to provide an initial impression of the simulation results and give a starting-point for model verification and validation.

Case Study: Predator-Prey - Basic Simulation and Simulation Results¶

Since the parameters used in the model are very abstract, it is difficult to select a meaningful initial parameter set and perform basic simulations. We estimate the average lifespan of a rodent to be approximately 1.5 years or 18 months, which results in a mortality rate of $\beta_1=1/18$ [deaths/month]. For fox we cannot use the same strategy since we cannot decouple their survival from hunting of prey. We further assume that female rodents have approximately two litters per year with 8 newborn, giving $\alpha_1=16/(12\dot 2)$ [births/month]. Without hunting, foxes cannot reproduce. Therefore, we assume $\alpha_2=0$ . Finally, we use $y_0=1000$ , as given by the data.

All other parameters, i.e. $\beta_2,\gamma,\delta$ and also $y_0$ , were initially guessed by iteratively running and visualising the model.

# fixed parameter values
alpha_1 = 8 / 12
beta_1 = 1 / 18
alpha_2 = 0.0
pred_0 = pred_data[0]
# free parameters - experiment with these!
beta_2 = 0.2
gamma = 0.002
delta = 0.1
prey_0 = pred_0 * 10
# initialise and run the model
sim = PredPreySim(alpha_1, beta_1, alpha_2, beta_2, gamma, delta)
t_model, prey_model, pred_model = sim.run(60, prey_0, pred_0)
# plot results
create_figure()
plot_prey_model(t_model, prey_model)
plot_pred_model(t_model, pred_model)
finish_plot()

Verification¶

In the next step, the model must be verified. This means that the implemented model is compared with the conceptual model to make sure that the model is correctly implemented. This process must be clearly distinguished from the validation process, in which the model or simulation results are compared with the real system and checks, whether the correct model was implemented. Typically, verification involves traditionally debugging up to the specification of test-scenarios for which a certain outcome is already deduced from formal model analysis.

Verification always has to come before validation (ex falso seuqitur quodlibet!).

Case Study: Predator-Prey - Verification¶

With the selected parameters, the simulation runs show oscillating behaviour, in which the rodent population always grows first. With a sufficiently large prey supply, the predator population then increases with a slight delay. As the number of predators grows, so does their food requirement, and the rodent population is pushed back. This then also affects the fox population, which also declines as the number of prey decreases. This allows the rodent population to slowly recover and the cycle begins anew. This behaviour was to be expected from the causal loop diagram, and, apparently, the numerical errors of to the used Runge-Kutta method do not disturb this.

In addition it is now possible to examine the behaviour of the model for $x_0=0$ and $y_0=0$ , which was investigated formally before. As expected, the simulation shows exponential growth for the prey and corresponding decline in the number of predators. We therefore consider the implemented model verified compared to the conceptual one. Apparently, the numerical errors due to the used Runge-Kutta method do not disturb the solution.

create_figure()
# run model without predators
t_model, prey_model, pred_model = sim.run(20, prey_0, 0)
plt.subplot(1, 2, 1)
# plot prey
plot_prey_model(t_model, prey_model)
finish_plot()
# run model without prey
plt.subplot(1, 2, 2)
t_model, prey_model, pred_model = sim.run(20, 0, pred_0)
# plot predators
plot_pred_model(t_model, pred_model)
finish_plot()

Validation and Validation Analysis¶

Validation is the process in which the model or model results are compared with reality. A clear distinction is made between qualitative and quantitative validation. While the former deals with the similarity of the model’s behaviour compared to the behaviour of the real system, the latter uses comparison metrics to compare the model’s outcomes with real data.

In both cases, a failed validation must be corrected. The first step is to evaluate whether the problem (a) lies with the model or the model structure, or (b) whether the parameterisation is insufficiently accurate.

Case Study: Predator-Prey - Validation and Validation Analysis¶

In the predator-prey example, a quantitative validation metric can be constructed using the given time series. The simulation can be evaluated for the given points in time and compared using, for example, a mean squared error.

def error(t_model, pred_model) -> float:
    """Computes the mean squared error between the given data and the predator time-series as given from a predator prey simulation
    :param t_model: reference time points
    :param pred_model: predator time series
    :return: mean squared error
    """
    inter = interp1d(t_model, pred_model)  # interpolate time series
    pred_model_inter = np.array([inter(t) for t in t_data])  # evaluate at data points
    diff = pred_model_inter - pred_data
    return sum(diff * diff) ** 0.5 / len(diff)


# set some free parameter values
beta_2 = 0.2
gamma = 0.001
delta = 0.1
prey_0 = 10000

# run the model
sim = PredPreySim(alpha_1, beta_1, alpha_2, beta_2, gamma, delta)
t_model, prey_model, pred_model = sim.run(60, prey_0, pred_0)

# plot the predators together with data
create_figure()
plot_pred_data()
plot_pred_model(t_model, pred_model)
plt.title(f"Error: {error(t_model, pred_model)}")
finish_plot()

While quantitative similarity is clearly given, still the model does not properly reflect the measurements. That means, at this point, the model cannot be considered valid and cannot yet be used for experimenting. Parameter values have to be refined first.

Parameter Identification¶

When validation fails due to invalid parameter values, they have to be adjusted. The corresponding refinement process is usually called parameter identification. There are basically two approaches here: (a) new data or information sources are collected in order to (re)determine the parameter values directly, or (b) parameter values are determined indirectly in such a way that the simulation result for certain reference scenarios matches the corresponding data. The latter process is also referred to as model calibration or model fitting. An existing validation metric can be used as an error function for the corresponding error-minimisation task. In case a large amount of training data is available (which is rather unusual for modelling and simulation applications) a train-test split, as is common in machine learning, can be performed to prevent overfitting.

Case Study: Predator-Prey - Parameter Identification¶

In the given setup, the four non-fixed parameters can be refined in a calibration setup using the reference data for the predators. This can be done either manually or by using automated optimisation.

In the following we simply state the indentified parameter values and do not go into any details on how they were found since this would go beyond the educational scope of this case study.

# identified parameter values
beta_2 = 0.24
gamma = 0.00092
delta = 0.13
prey_0 = 9900

# run the simulation and check fit
sim = PredPreySim(alpha_1, beta_1, alpha_2, beta_2, gamma, delta)
t_model, prey_model, pred_model = sim.run(60, prey_0, pred_0)

# plot predator simulation together with data
create_figure()
plot_pred_data()
plot_pred_model(t_model, pred_model)
plt.title(
    f"Params: {beta_2=:.02f},{gamma=:.05f},{delta=:.02f},{prey_0=:.02f}\nError: {error(t_model, pred_model)}"
)
finish_plot()

The fit has clearly improved and would now be acceptable for experiements. However, the model unfortunately fails replicating the observed behaviour qualitatively. Due to its nature, the model only provides perfectly periodic solutions, but the data show a temporary drop in the peaks.

The validation failed again. However, this time it cannot be fixed by a simple parameter identification process. The conceptual model structure is apparently not ok and must be redesigned.

Case-Study: Predator-Prey (Re-Modelling)¶

The modeller is now searching for system components that are not accurately reflected in the defined conceptual model and could be responsible for the observed qualitative differences. He recognises a possible weakness in the model in the fact that prey animals grow unbounded in the absence of predators, which is not realistic.

In order to model this behaviour more realistically, the modeller introduces a resource cap for the growth of prey animals. In the causal loop diagram or stock and flow diagram, this change is shown as depicted in the figure below (blue additions).

In the differential equation model, this modification leads to the following equation for the prey animals:

\frac{dx}{dt}(t)=x(t)\left(\alpha_1\cdot(1-\xi(t))\frac{M-x(t)}{M}-\beta_1\right)-x(t)y(t)\gamma.

(5)

Hereby, $M$ [prey] stands for the resource limit of the prey animals e.g. due to internal competition.

In the absence of predators, i.e. $y(t)\equiv 0$ , and without regarding the intervention, the equation states

\frac{dx}{dt}(t)=x(t)\left(\alpha_1\frac{M-x(t)}{M}-\beta_1\right) = x(t)\left((\alpha_1-\beta_1)-x(t)\frac{\alpha_1}{M}\right).

(6)

As a result, the growth is no longer exponential, but logistic, approaching an upper bound when $t\rightarrow \infty$ .

We will furthermore call this model damped predator-prey model.

Case-Study: Predator-Prey (Re-Programming)¶

In the next step, we adapt the implementation. To make sure, that only the model equations change, we simply extend the class implemented earlier.

class PredPreyDampedSim(PredPreySim):
    def __init__(
        self,
        alpha_1: float,
        beta_1: float,
        alpha_2: float,
        beta_2: float,
        gamma: float,
        delta: float,
        M: float,
        xi: float = 0.0,
        xi_time: [float, float] = None,
    ):
        """Initialises a new damped predator-prey simulation model
        :param alpha_1: birth rate prey
        :param beta_1: death rate prey
        :param alpha_2: birth rate predators
        :param beta_2: death rate predators
        :param gamma: number of prey anmials eaten per prey and per predator animal
        :param delta: number of newborn predator from one successful prey hunt
        :param M: resource limit for the prey animals
        :param xi: effectiveness of the chemical intervention on the prey's birth rate
        :param xi_time: timeframe in which the intervention is active
        """
        super().__init__(alpha_1, beta_1, alpha_2, beta_2, gamma, delta, xi, xi_time)
        self.M = M  # the model requires one additional parameter

    def rhs(self, t: float, xy: np.ndarray) -> np.ndarray:
        """Right hand side of the differential equation.
        Overrides the one from the undamped model.
        :param t: current time
        :param xy: current number of prey and predator animals as vector
        :return: vector of change of the prey and predator animals
        """
        prey = xy[0]
        pred = xy[1]
        prey_damping = (self.M - prey) / self.M  # this is the new term!
        eaten = pred * prey * self.gamma
        newpred = eaten * self.delta
        if self.xi_time[0] < t <= self.xi_time[1]:
            preyBirth = (
                prey * self.alpha_1 * (1 - self.xi) * prey_damping
            )  # multiply with new factor
        else:
            preyBirth = prey * self.alpha_1 * prey_damping  # multiply with new factor
        preyDeath = prey * self.beta_1 + eaten
        predBirth = pred * self.alpha_2 + newpred
        predDeath = pred * self.beta_2
        return np.array([preyBirth - preyDeath, predBirth - predDeath])

Case-Study: Predator-Prey (Re-Validation and Validation Analysis)¶

We perform several new validation experiments with the new model to assess, whether it is now capable of producing the expected behaviour, i.e. that the peaks of the predator animals decline with time.

# some experimental parameter values
beta_2 = 0.2
gamma = 0.001
delta = 0.1
prey_0 = pred_0 * 10
M = pred_0 * 50

# run damped and old model and compare the results
sim = PredPreyDampedSim(alpha_1, beta_1, alpha_2, beta_2, gamma, delta, M)
sim_old = PredPreySim(alpha_1, beta_1, alpha_2, beta_2, gamma, delta)
t_model, prey_model, pred_model = sim.run(200, prey_0, pred_0)
t_model_old, prey_model_old, pred_model_old = sim_old.run(200, prey_0, pred_0)

# make plots
create_figure()
plot_prey_model(t_model, prey_model, label="rodents (damped simulation)")
plot_prey_model(
    t_model_old,
    prey_model_old,
    label="rodents (undamped simulation)",
    linestyle="dashed",
)
plot_pred_model(t_model, pred_model, label="foxes (damped simulation)")
plot_pred_model(
    t_model_old, pred_model_old, label="foxes (undamped simulation)", linestyle="dashed"
)
finish_plot()

Case-Study: Predator-Prey (Re-Identify Parameters)¶

Indeed the damped model is now capable of producing the expected qualitative behaviour. Thus we may perform an new parameter identification round to fit the model to the data.

Again, we do not go into any details on how these values have been found but simply show the results of the identification process.

# identified parameter values
beta_2 = 0.23
gamma = 0.00066
delta = 0.06
prey_0 = 35000
M = 170000

# run the simulation and check fit
sim = PredPreyDampedSim(alpha_1, beta_1, alpha_2, beta_2, gamma, delta, M)
t_model, prey_model, pred_model = sim.run(60, prey_0, pred_0)

# plot predator results together with data
create_figure()
plot_pred_data()
plot_pred_model(t_model, pred_model)
plt.title(
    f"Params: {beta_2=:.02f},{gamma=:.05f},{delta=:.03f},{prey_0=:.0f},{M=:.0f}\nError: {error(t_model, pred_model)}"
)
finish_plot()

Experiments with Model¶

Once the model has passed the validation process, it is ready for the actual simulation studies. For the modeller, this means that they can finally reap the rewards of all the hard work they have put into creating and validating the model. However, it is important to note that the actual model experiments must be carried out with sufficient care, as there are a number of key points to consider.

Firstly, it is crucial that the experiments are conducted in a reproducible and comprehensible manner. Although this applies to the entire process described so far, the model results on which the final results of the entire study are based must be documented particularly accurately and in detail. A properly specified and documented experimental setup, which allows to precisely replicate the results, is the cornerstone of this documentation process.

Furthermore, it is essential to consider the uncertainties inherent in the system itself on the one hand, and those associated with the modelling process and its simplifications on the other. These should also be quantified, e.g. via variation of parameters, inputs or model-inherent stochasticity.

Finally, the results must also be presented in such a way that they are comprehensible and understandable to decision-makers. Only then will they be useful to them.

Case Study: Predator-Prey - Experiments with Model¶

With the parameters identified for the system without intervention, intervention scenarios can now be calculated. To do this, the model is extrapolated from the last 60 months into the future for a further 120 months. At fixed points in time, a chemical pesticide is then applied, which reduces the birth rate of rodents by 90% for two months at a time, i.e. $\xi=0.9$ .

In the first experiment, we compare the model results for four different points in time for the intervention: months $80,90,100$ and 110.

create_figure()
xi = 0.9  # effectiveness of the intervention
di = 2  # duration of the intervention
ind = 0
for ti in [80, 90, 100, 110]:  # vary the intervention month
    ind += 1
    plt.subplot(2, 2, ind)
    # initialise simulation with specified intervention month
    sim = PredPreyDampedSim(
        alpha_1, beta_1, alpha_2, beta_2, gamma, delta, M, xi, [ti, ti + di]
    )
    t_model, prey_model, pred_model = sim.run(180, prey_0, pred_0)
    plot_prey_model(t_model, prey_model, zorder=1)
    # highlight fitting and intervention period
    yl = plt.ylim()
    plt.fill_between(
        [0, 60],
        [0, 0],
        [yl[1], yl[1]],
        color="k",
        alpha=0.3,
        label="historic",
        zorder=0,
    )
    plt.fill_between(
        [ti, ti + di],
        [0, 0],
        [yl[1], yl[1]],
        color="r",
        alpha=0.3,
        label="intervention period",
        zorder=0,
    )
    plt.ylim(0, yl[1])
    plt.xlim(0, 180)
    plt.title(f"intervention at month {ti}")
    finish_plot()
plt.tight_layout()

The results are highly interesting to interpret. When the intervention is activated in a period with a high number of rodents, then the rodent population is damped and becomes more stable afterwards. If, however, the intervention is activated in a period in which the number is low, it will lead to even more fluctuating population figures with even higher population peaks.

The question remains as to how reliable statements about the timing of the intervention are, e.g. month 80 looks very promising now. To assess this, uncertainty in parameterisation and input can be investigated using random variables.

In the second experiment we vary the predator death rate $\beta_2$ and the initial prey count $x_0$ by 10% and activate the intervention at month 80.

ti = 80  # intervention month
di = 2
vary = 0.1  # uncertainty range (10%)
create_figure()
cmp = plt.get_cmap("jet")
for i in range(100):
    f1 = np.random.random() * 2 - 1.0  # produces a random number between -1 and 1
    f2 = np.random.random() * 2 - 1.0  # produces a random number between -1 and 1
    sim = PredPreyDampedSim(
        alpha_1,
        beta_1,
        alpha_2,
        beta_2 * (1 + vary * f1),
        gamma,
        delta,
        M,
        xi,
        [ti, ti + di],
    )
    t_model, prey_model, pred_model = sim.run(180, prey_0 * (1 + vary * f2), pred_0)
    if i == 0:  # only need label once
        lbl = "rodents (simulation)"
    else:
        lbl = None
    plot_prey_model(
        t_model, prey_model, zorder=1, linewidth=0.3, color=cmp(i / 100), label=lbl
    )
# highlight fitting and intervention timespan
yl = plt.ylim()
plt.fill_between(
    [0, 60], [0, 0], [yl[1], yl[1]], color="k", alpha=0.3, label="historic"
)
plt.fill_between(
    [ti, ti + di],
    [0, 0],
    [yl[1], yl[1]],
    color="r",
    alpha=0.3,
    label="intervention period",
)
plt.ylim(0, yl[1])
plt.xlim(0, 180)
plt.title(f"intervention at month {ti}")
finish_plot()

Problem Solution¶

Finally, it remains to solve the original problem using the new insights gained from the simulation experiments. It should be noted that this must be done using an “inverse abstraction process” that takes into account how the respective simulation results are to be interpreted in the real system (a model is not the reality!). Note that this process must be carried out by the modeller, not by the programmer or the decision maker, because only they are familiar with the abstraction process involved in model development.

In this final step, one must also be prepared for the fact that often some of the original questions can not be answered or must be answered differently than originally expected due to the limitations of the model. At the same time, the experiments may provide answers to questions that were not even asked at the beginning. This requirement for openness to results makes modelling and simulation studies very difficult to plan from an exonomic perspective on the one hand, but also very exciting and enriching on the other.

Case-Study: Predator-Prey - Problem Solution¶

The first experiment clearly showed that the success of the intervention depends heavily on timing. If the wrong time is chosen, the intervention can disrupt the balance between predator and prey and ultimately cause even more damage than without it, due to higher population peaks.

However, the second experiment shows that the simulation model cannot help to find a good intervention time, as the inherent uncertainty of the system is too great. Even a deviation of 10% in parameter values can mean that a time initially identified as favourable is actually unfavourable. However, the model has shown that deployment is useful when the number of rodents is as high as possible.

Finally, experiments with the model have shown that a predator-prey dynamics as the one observed between the invasive rodents and the foxes will generally tend towards a constant equilibrium behaviour with negligibly small population peaks - that is, if it is not continuously disturbed by external influences not yet regarded in the model, e.g. climate, weather, soil, etc.. A short-term intervention in the system, however, such as the planned use of chemicals cannot change the dynamics between the two species in the long term, but can at most delay or accelerate them. This should be taken into account when considering the one-off use of chemicals.

References¶

Banks, J. (2005). Discrete event system simulation. Pearson Education India.
Sargent, R. G. (2010). Verification and validation of simulation models. Proceedings of the 2010 Winter Simulation Conference, 166–183. 10.1109/wsc.2010.5679166