Introduction to the python library for system dynamics, that is part of the BPTK-Py business simulation framework.

# SD DSL: Under The Hood

Creating A Domain Specific Language For System Dynamics Simulations

Part of the BPTK Framework is the System Dynamics Domain Specific Language (SD DSL). Using the SD DSL, you can easily build System Dynamics Mmodels directly in Python using an intuitive Syntax.

This notebook takes a look under the hood of the SD DSL and explains how it creates the underlying equations that are needed to run a simulation.

This should give you a better understanding of System Dynamics, of the SD DSL and of Python metaprogramming techniques.

The following diagram contains the three key System Dynamics language elements:

• Stocks
• Flows
• Converters

At any given time $$t$$, a stock is equal to the value of the stock at time $$t-1$$ plus the sum of all inflows, minus the sum of all outflow at time $$t-1$$ are simply the sum of the flows that flow into and out of them, over time.

Flows and converters are simply functions of their inputs.

The general mathematical equations for the stocks above are:

$\begin{equation*} stock(t)=stock(t-1)+\sum_{inflows}inflow(t-1)-\sum_{outflows}outflow(t-1) \end{equation*}$

$\begin{equation*} flow(t) = function(input_1,...,input_n) \end{equation*}$

$\begin{equation*} converter(t) = function(input_1,...,input_n) \end{equation*}$

The concrete equations for the model in the diagram above are:

$\begin{equation*} stock(t)=stock(t-1)+\sum_{inflows}inflow(t-1) \end{equation*}$

$\begin{equation*} flow(t) = stock(t)*rate(t) \end{equation*}$

$\begin{equation*} rate(t) = 0.1 \end{equation*}$

$\begin{equation*} stock(0) = 1 \end{equation*}$

The SD DSL allows us to write SD models in Python using a simple syntax:

stock = model.stock("stock")
flow = model.flow("flow")
rate = model.converter("rate")
stock.initial_value = 1
rate.equation=0.1
flow.equation=stock*rate
stock.equation=flow

But under the hood, the SD DSL needs to build the equations listed above to ensure that we can simulate the model.

Let’s take a look at how this is done, step by step.

## Step 1: Write Python Code For The Equations Above

The SD DSL uses lambda functions lambda functions to implement the equation if the System Dynamics model.

import pandas as pd
import matplotlib
stock = lambda t: stock(t-1) + flow(t-1) if t>0 else 1
flow = lambda t: stock(t) * rate(t)
rate = lambda t: 0.1
stock(10)
2.5937424601

## Step 2: Plot The Resulting Equation

Plot the stock equation for 20 timesteps, e.g. using a pandas dataframe.

df = pd.DataFrame((stock(t)) for t in range(0,20) )
df.plot(kind="line")
<Axes: >

## Step 3: Keep Equations in A Dictionary

Assume we want to keep equations in a “model” and add new ones dynamically: create a dictioanry of equations that are accessed by equation names as strings, e.g. equation["stock"](5)

equations={}
equations["stock"] = lambda t: equations["stock"](t-1) + equations["flow"](t-1) if t>0 else 1
equations["flow"] = lambda t: equations["stock"](t) * equations["rate"](t)
equations["rate"] = lambda t: 0.1
equations["stock"](10)
2.5937424601

## Step 4: Use Memoization To Improve Performance

In order to calculate the value of a stock at a given time $$t$$, you need to know all preceding values. The way the model is encoded at the moment, this means that we first calculate the value at time $$t = 1$$, then to calculate the value at time $$t = 2$$ we need to calcualte the values at time $$t=1$$ and then at time $$t =2$$ and so on … this means that even models with only few timesteps can take a very long time to evaluate.

The way around this of course is to remember the value of each equation at each time step, so that you don’t have to recurse through all timesteps at every timestep.

We do this using a technique called memoization.

%timeit equations["stock"](10)
151 µs ± 1.08 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
%timeit equations["stock"](20)
519 ms ± 3.33 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
equations={}
memo={}

The memoization function stores the value of each equation at each timestep. When asking for the value of an equation at a given time, it first checks if the value has already been calculated. If yes, it returns that value. If no, it calculates the value, stores it and then returns it.

def memoize(equation, time):

try:
mymemo = memo[equation]
except:
# In case the equation does not exist in memo
memo[equation] = {}
mymemo = memo[equation]
if time in mymemo.keys():
return mymemo[time]
else:
result = equations[equation](time)
mymemo[time] = result

return result
equations["stock"] = lambda t: memoize("stock",t-1) + memoize("flow",t-1) if t>0 else 1
equations["flow"] = lambda t: memoize("stock",t) * memoize("rate",t)
equations["rate"] = lambda t: 0.1
equations["stock"](5)
1.61051
%timeit equations["stock"](20)
187 ns ± 0.171 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
%timeit equations["stock"](30)
187 ns ± 1.23 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
equations["stock"](50)
117.39085287969533
df = pd.DataFrame((equations["stock"](t)) for t in range(0,50) )
df.plot(kind="line")
<Axes: >

## Step 5: Compile Equations From Strings

Now that we know what equations we want, we need a find a way of generating them directly. We do this by first generating the code in string form and then compiling that string into code using the evalfunction.

equations["stock"]=eval("lambda t: memoize('stock',t-1) + memoize('flow',t-1) if t>0 else 1")
equations["stock"](5)
1.61051
equations["stock"]
<function __main__.<lambda>(t)>

## Step 6: Encapsulate The Basic Functionality In Classes

Now that we have all the building blocks for our SD DSL, we would like to encapsulate them in classes that allow us to build a model at runtime.

We want to be able to build a model using objects at runtime, the underlying equations should be generated for us.

• Create classes for the model (that will hold the model elements and equations) and for the elements themselves.
• New elements are created using the model as a factory (e.g. another_stock = model.stock("name of another stock"))
• Elements know which model they belong to
• The equations for each element are defined using an equation method on the elements (e.g. stock.equation=flow)
• The memoization function should be part of the model class
• Each element should have a term function element.term that generate a string representation of itself (e.g model.memoize("stock",t))
• Remember that stocks have an initial value
class Model:
def __init__(self):
self._stocks = {}
self._flows = {}
self._converters = {}
self._memo={}

def stock(self, name):
if name not in self._stocks:
self._stocks[name] = Stock(model, name)
return self._stocks[name]

def flow(self, name):
if name not in self._flows:
self._flows[name] = Flow(model, name)
return self._flows[name]

def converter(self, name):
if name not in self._converters:
self._converters[name] = Converter(model, name)
return self._converters[name]

def memoize(self,equation, arg):

try:
mymemo = self._memo[equation]
except:
# In case the equation does not exist in memo
self._memo[equation] = {}
mymemo = self._memo[equation]
if arg in mymemo.keys():
return mymemo[arg]
else:
result = self._equations[equation](arg)
mymemo[arg] = result

return result
class Element:
def __init__(self, model, name):
self._model = model
self._name = name
self._equation = None

@property
def equation(self):
return self._equation

@equation.setter
def equation(self, equation):
self._equation = equation

def term(self,time="t"):
return "model.memoize('{}',{})".format(self._name, time)
class Stock(Element):

def __init__(self,model, name):
super(Stock, self).__init__(model, name)
self._initial_value = 0

@property
def initial_value(self):
return self._initial_value

@initial_value.setter
def initial_value(self, value):
self._initial_value=value

class Flow(Element):
pass

class Converter(Element):
pass

model = Model()
stock = model.stock("stock")
stock.term()
"model.memoize('stock',t)"
stock.initial_value=1
rate = model.converter("rate")
rate.term()
"model.memoize('rate',t)"
flow = model.flow("flow")
flow.term()
"model.memoize('flow',t)"
stock.equation = flow
stock._equation
<__main__.Flow at 0x12e496b90>
rate.equation=10

## Step 7: Auto-Generate The Equations

Extend the framework to auto-generate the underlying lambda-equations and add them to the model. For the first version assume we just have a stock and a constant flow.

Override the call operator to evalue the element at a particular timestep (e.g. stock(10))

Create a plot function that evaluates an element in the range 0 to 100.

class Model:
def __init__(self):
self._stocks = {}
self._flows = {}
self._converters = {}
self._memo={}
self._equations={}

def stock(self, name):
if name not in self._stocks:
self._stocks[name] = Stock(model, name)
return self._stocks[name]

def flow(self, name):
if name not in self._flows:
self._flows[name] = Flow(model, name)
return self._flows[name]

def converter(self, name):
if name not in self._converters:
self._converters[name] = Converter(model, name)
return self._converters[name]

def memoize(self,equation, arg):

try:
mymemo = self._memo[equation]
except:
# In case the equation does not exist in memo
self._memo[equation] = {}
mymemo = self._memo[equation]
if arg in mymemo.keys():
return mymemo[arg]
else:
result = self._equations[equation](arg)
mymemo[arg] = result

return result
class Element:
def __init__(self, model, name):
self._model = model
self._name = name
self._equation = None
self._function_string=""

@property
def equation(self):
return self._equation

@equation.setter
def equation(self, equation):
self._equation = equation
self._function_string = "lambda model, t: {}".format(equation)
self.generate_function()

def generate_function(self):
"""
Generate the function using the function_string value and eval()
:return: None
"""
fn = eval(self._function_string)
self._model._equations[self._name] = lambda t: fn(self._model, t)
self._model._memo[self._name] = {}

def term(self,time="t"):
return "model.memoize('{}',{})".format(self._name, time)

def plot(self):
df = pd.DataFrame((self(t)) for t in range(0,100) )
df.plot(kind="line")

def __call__(self, *args, **kwargs):
return self._model._equations[self._name](args[0])

def __str__(self, other):
return self.term()
class Stock(Element):
type="Stock"

def __init__(self,model, name):
super(Stock, self).__init__(model, name)
self._initial_value = 0

@property
def initial_value(self):
return self._initial_value

@initial_value.setter
def initial_value(self, value):
self._initial_value=value

@property
def equation(self):
return self._equation

@equation.setter
def equation(self, equation):
self._equation = equation

self._function_string = "lambda model, t : ( ("
self._function_string += str(self._initial_value)
self._function_string += ") if (t <= 0) else (model.memoize('{}',t-1)) ".format(self._name)

if self._equation is not None:
self._function_string += "+ ("
self._function_string += self._equation.term("t-1")
self._function_string += ") )"
else:
self._function_string += ")"

self.generate_function()

class Flow(Element):
type="Flow"

class Converter(Element):
type="Converter"

model = Model()
stock = model.stock("stock")
stock.term()
"model.memoize('stock',t)"
stock.initial_value=1
flow = model.flow("flow")
flow.term()
"model.memoize('flow',t)"
flow.equation=1
stock.equation = flow
stock._equation
<__main__.Flow at 0x10ee78748>
model._equations
{'flow': <function __main__.Element.generate_function.<locals>.<lambda>(t)>,
'stock': <function __main__.Element.generate_function.<locals>.<lambda>(t)>}
flow(10)
1
stock(10)
11
stock.plot()

## Step 8: Extend the Equation Syntax Using Operators

In order to re-create our initial model, we need to formulate an equation such as flow.equation = stock*rate.

Override the __mul__ methods of the Element and create a multiplication operator to deal with this.

class Model:
def __init__(self):
self._stocks = {}
self._flows = {}
self._converters = {}
self._memo={}
self._equations={}

def stock(self, name):
if name not in self._stocks:
self._stocks[name] = Stock(model, name)
return self._stocks[name]

def flow(self, name):
if name not in self._flows:
self._flows[name] = Flow(model, name)
return self._flows[name]

def converter(self, name):
if name not in self._converters:
self._converters[name] = Converter(model, name)
return self._converters[name]

def memoize(self,equation, arg):

try:
mymemo = self._memo[equation]
except:
# In case the equation does not exist in memo
self._memo[equation] = {}
mymemo = self._memo[equation]
if arg in mymemo.keys():
return mymemo[arg]
else:
result = self._equations[equation](arg)
mymemo[arg] = result

return result
class Element:
def __init__(self, model, name):
self._model = model
self._name = name
self._equation = None
self._function_string=""

@property
def equation(self):
return self._equation

@equation.setter
def equation(self, equation):
self._equation = equation
self._function_string = "lambda model, t: {}".format(equation)
self.generate_function()

def generate_function(self):
"""
Generate the function using the function_string value and eval()
:return: None
"""
fn = eval(self._function_string)
self._model._equations[self._name] = lambda t: fn(self._model, t)
self._model._memo[self._name] = {}

def term(self,time="t"):
return "model.memoize('{}',{})".format(self._name, time)

def plot(self):
df = pd.DataFrame((self(t)) for t in range(0,100) )
df.plot(kind="line")

def __call__(self, *args, **kwargs):
return self._model._equations[self._name](args[0])

def __mul__(self, other):
return MultiplicationOperator(self, other)

def __str__(self):
return self.term()
class MultiplicationOperator():

def __init__(self, element_1, element_2):
self._element_1 = element_1
self._element_2 = element_2

def term(self, time="t"):
return "(" + self._element_1.term(time) + ") * (" + self._element_2.term(time) + ")"

def __str__(self):
return self.term()
class Stock(Element):
type="Stock"

def __init__(self,model, name):
super(Stock, self).__init__(model, name)
self._initial_value = 0

@property
def initial_value(self):
return self._initial_value

@initial_value.setter
def initial_value(self, value):
self._initial_value=value

@property
def equation(self):
return self._equation

@equation.setter
def equation(self, equation):
self._equation = equation

self._function_string = "lambda model, t : ( ("
self._function_string += str(self._initial_value)
self._function_string += ") if (t <= 0) else (model.memoize('{}',t-1)) ".format(self._name)

if self._equation is not None:
self._function_string += "+ ("
self._function_string += self._equation.term("t-1")
self._function_string += ") )"
else:
self._function_string += ")"

self.generate_function()

class Flow(Element):
type="Flow"

class Converter(Element):
type="Converter"

model = Model()
stock = model.stock("stock")
flow = model.flow("flow")
rate = model.converter("rate")
rate.equation=0.1
stock.initial_value=1
stock.equation = flow
flow.equation = stock*rate
stock(3)
1.3310000000000002
stock.plot()