A Simple Python Library For System Dynamics

Building System Dynamics simulations interactively in Jupyter using Python

Oliver Grasl
Monday, October 1, 2018

An overview of how to build a System Dynamics model using the BPTK-Py framework

Build System Dynamics Simulations Interactively in Jupyter Using Python

We love building computational models and our favorite environments for this kind of explorative, analytical work are Jupyter and Python.

To make computational modeling easier we are developing a simple framework that currently supports System Dynamics and Agent-based modeling.

We first introduced the BPTK-Py framework in our recent post Writing Computational Essays Based on Simulation Models.

Since then we have created some new functionality that allows you to build System Dynamics models and Agent-based models interactively in Jupyter using Python. To make the model building as simple as possible, we have created a simple, domain-specific language (DSL) that supports both System Dynamics and Agent-based modeling and hides much of the underlying complexity of computational models.

This language not only allows you to create System Dynamics models and Agent-based models, but you can also even mix the two to create “hybrid” simulation models.

Having such a DSL is useful for several reasons:

Build models interactively in Jupyter, making the modeling process very effective.
Python novices can focus on the modeling, without needing to know much about Python.
Python experts can mix their models with other analytical frameworks, e.g. Machine Learning toolkits.

Needless to say, the new functionality seamlessly integrates with the rest of the BPTK-Py framework, so you can use all the high-level scenario management and plotting functions that are part of the framework.

In this post I focus on how to build a System Dynamics model using the framework, I will take a look at Agent-based modeling in a future post.

The post is also available as a Jupyter notebook in our BPTK-Py Tutorial, which is available for download here.

A Simple Model to Demonstrate The Library

To illustrate the DSL, we will build the simple project management model we introduced in our Introduction to System Dynamics.

The project management model is really simple and just contains a few stocks, flows and converters, as you can see in the following diagram:

So even if you do not know the model you should be able to follow this post very easily.

To get started, we first need to import the library and in particular the SD function library into our notebook.

from BPTK_Py import Model
from BPTK_Py import sd_functions as sd

The SD function library contains the functions and operators needed to define model equations (these are the built-ins you will know from your visual modeling environments, such as Stella or Vensim). Because the library contains functions such as

min

and

max

, which are also part of the Python standard library, we import the SD function library with the prefix sd to avoid naming conflicts.

Next, we create a model using the

Model

class. Our model will contain all our model elements such as stocks, flows, converters and constants.

model = Model(starttime=0,stoptime=120,dt=1,name='SimpleProjectManagament')

Creating model elements is really easy:

openTasks = model.stock("openTasks")

As you can see our convention is to use the camel casing naming convention for the model elements and to create Python variables for the elements that carry the same name.

Once a model element has been defined in this manner, we only need to refer to the Python variable and do not need to reference the model element (i.e. we can use

openTasks

in our equations, as opposed to using

model.stock("openTasks")

). This saves a lot of typing.

Let us define the other model elements and variables now too so that we can then concentrate on the equations:

closedTasks = model.stock("closedTasks")
staff = model.stock("staff")
completionRate = model.flow("completionRate")
currentTime = model.converter("currentTime")
remainingTime = model.converter("remainingTime")
schedulePressure = model.converter("schedulePressure")
productivity = model.converter("productivity")
deadline = model.constant("deadline")
effortPerTask = model.constant("effortPerTask")
initialStaff = model.constant("initialStaff")
initialOpenTasks = model.constant("initialOpenTasks")

Note that in our models we differentiate between constants and converters – this is not strictly necessary from a System Dynamics point of view, but it makes it easier to check the model for errors.

Now let us initialize our stocks – to do this, we just need to set the

initial_value

property of the stocks. The initial value can either be a numerical constant or a constant element.

closedTasks.initial_value = 0
staff.initial_value = initialStaff
openTasks.initial_value = initialOpenTasks

Defining the model equations is really easy: each model variable has an equation property, the equation itself is written much as you would in a visual modeling environment, using the other model variables as necessary.

Defining constants is particularly easy:

deadline.equation = 100
effortPerTask.equation = 1
initialStaff.equation = 1
initialOpenTasks.equation = 100

The

currentTime

variable tracks the simulation time, which is captured by the

time

function in the SD function library.

currentTime.equation=sd.time()

The

remainingTime

equals the difference between the

deadline

and the

currentTime

remainingTime.equation = deadline - currentTime

So you see, thanks to the DSL defining equations is very intuitive!

The equations for the stocks are also really simple – they just contain the inflows (with a positive sign) and the outflows (with a negative sign).

A quick look at the diagram above tells us that the

openTasks

only have an outflow (defined by the

completionRate

) and the

closedTasks

only have an inflow (also defined by the completionRate):

openTasks.equation = -completionRate
closedTasks.equation = completionRate

The

schedulePressure

is a dimensionless ratio that compares the required effort to complete all remaining open tasks to the remaining work capacity.

We use the

min

and

max

functions from the SD function library to ensure that no division by zero occurs and that the schedule pressure is capped at 2.5:

schedulePressure.equation = sd.min(
   (openTasks*effortPerTask)/
   (staff*sd.max(remainingTime,1)),
   2.5
)

We define the productivity in our model using a non-linear relationship (depending on the schedule pressure). We capture this relationship in a lookup table that we store in the

points

property of the model (using a Python list):

model.points["productivity"] = [
    [0,0.4],
    [0.25,0.444],
    [0.5,0.506],
    [0.75,0.594],
    [1,1],
    [1.25,1.119],
    [1.5,1.1625],
    [1.75,1.2125],
    [2,1.2375],
    [2.25,1.245],
    [2.5,1.25]
]

We can easily plot the lookup table to see whether it has the right shape:

model.plot_lookup("productivity")

A Simple Python Library for System Dynamics graph 1

The productivity equation is then defined via a lookup function – in our case

productivity

depends non-linearly on

schedulePressure

as defined in the lookup table:

productivity.equation = sd.lookup(schedulePressure,"productivity")

The last equation we need to define is that of the

completionRate

– the completion rate is defined by the number of staff working on the project divided by the effort per task. We then multiply this by the (average) productivity of the staff. The completion rate may never be larger than the number of

openTasks

, so we constrain it using the

min

function.

completionRate.equation = sd.max(0,
   sd.min(openTasks, 
      staff*(productivity/effortPerTask)
   )
)

Now that we have defined all necessary equations, we are ready to run the model. The easiest way is to evaluate a model variable at a particular timestep – this approach is particularly useful if you are building the model interactively (e.g. in a Jupyter notebook) and you want to test intermediate results.

closedTasks(80), closedTasks(100), closedTasks(120)

(80.0, 100.0, 100.0)

Let us play with different settings for the deadline:

deadline.equation = 120
closedTasks(80), closedTasks(100), closedTasks(120)

(63.33020661244643, 81.06644489208418, 99.99777243819346)

deadline.equation=80
closedTasks(80), closedTasks(100), closedTasks(120)

(92.6853060260874, 100.00000000000004, 100.00000000000004)

Of course, we can also plot the variables in a graph straight away using the element’s

plot()

method.

closedTasks.plot()

A Simple Python Library for System Dynamics graph 2

Now that we have a model, we can use the powerful scenario management built into the BPTK-Py framework.

To do that, we first need to instantiate the framework:

import BPTK_Py
bptk = BPTK_Py.bptk()

Then we set up a scenario manager using a Python dictionary. The scenario manager identifies the baseline constants of the model:

scenario_manager = {
    "smSimpleProjectManagementDSL":{
    
    "model": model,
    "base_constants": {
        "deadline": 100,
        "initialStaff": 1,
        "effortPerTask": 1,
        "initialOpenTasks": 100,

    },
    "base_points":{
            "productivity": [
                [0,0.4],
                [0.25,0.444],
                [0.5,0.506],
                [0.75,0.594],
                [1,1],
                [1.25,1.119],
                [1.5,1.1625],
                [1.75,1.2125],
                [2,1.2375],
                [2.25,1.245],
                [2.5,1.25]
            ]
    }
 }
}

The scenario manager has to be registered as follows:

bptk.register_scenario_manager(scenario_manager)

Once we have this, we can define and register (one or more) scenarios as follows:

bptk.register_scenarios(
    scenarios =
        {
            "scenario80": {
                "constants": {
                    "initialOpenTasks": 80
                }
            }
        }
    ,
    scenario_manager="smSimpleProjectManagementDSL")

We can then plot the scenario as follows:

bptk.plot_scenarios(
    scenarios="scenario80",
    scenario_managers="smSimpleProjectManagementDSL",
    equations="openTasks")

A Simple Python Library for System Dynamics graph 3

Let us register a few more scenarios:

bptk.register_scenarios(
    scenarios =
    {
         "scenario100": {
        
        },
        "scenario120": {
            "constants": {
                "initialOpenTasks" : 120
            }
        }
    },
    scenario_manager="smSimpleProjectManagementDSL")

scenario100

is equivalent to the base settings, hence we do not need to define any settings for it.

Now we can easily compare the scenarios:

bptk.plot_scenarios(
    scenarios="scenario80,scenario100,scenario120",
    scenario_managers="smSimpleProjectManagementDSL",
    equations="openTasks",
    series_names={
        "smSimpleProjectManagementDSL_scenario80_openTasks":"scenario80",
        "smSimpleProjectManagementDSL_scenario100_openTasks":"scenario100",
        "smSimpleProjectManagementDSL_scenario120_openTasks":"scenario120"
    }
)

A Simple Python Library for System Dynamics graph 4

This completes our quick tour of the SD DSL within the Business Prototyping Toolkit. The BPTK Framework is available under the MIT License on PyPi, so you can start using it right away.

You can download a tutorial (which includes this blog post as a Jupyter notebook) on the Business Prototyping Toolkit.

The tutorial contains also illustrates some more advanced techniques, in particular also on how you can use the SD DSL in Python without using Jupyter.

Conclusion

This post introduced a simple domain-specific language for System Dynamics, implemented in Python. It lets you create System Dynamics in Python and supports interactive modeling in Jupyter.

Creating System Dynamics models directly in Python is particularly useful if you have the need to extend your SD models with your own SD functions or you would like to combine your models with other computational models such as Agent-based models or mathematical models.

Writing Computational Essays Based on Simulation Models Overview Using Simulations to Train And Test Machine Learning Applications