In [1]:
import pandas


The system dynamics of diapers with weekly pickup/delivery::

  q  +-----+  r  +-----+  s
---->|  C  |---->|  D  |---->
^^  +-----+     +-++--+
||                ||
++= = = = = = = = ++


$C$ = stock of clean diapers

$D$ = stock of dirty diapers

$q$ = in-flow of clean diapers

$r$ = flow of clean diapers to dirty diapers

$s$ = out-flow of dirty diapers

To make this totally precise we need a system of equations relating the stocks and flows. In this case, a system of difference equations with a timestep of one day works fine:

$C(t+1) = C(t) + q(t) - r(t)$

$D(t+1) = D(t) + r(t) - s(t)$

So far this seems like an example from intro applied math, and 100 other subjects. But what makes it truely a system dyamics modeling example is the dashed arrow from compartment $D$ to flow $q$. This signifies feedback, a core feature of SDM. It is actually what makes me want this simulation, because it makes things complicated (especially since I am a new parent, and therefore I am not sleeping much!) The in-flow of clean diapers happens weekly (on Fridays... don't forget to take the dirty ones out on Thursday!), and the number of diapers that flow in is a function of the number of diapers that flow out... a week ago. In equations, $$q(t) = s(t-7) + A(t) - C(t-7).$$

I've added something new here, $A(t)$ which is the stock of diapers I ask for. It started out at 70, but I was worried about running out, so I raised it to 100.

So what is left to specify? $r(t)$ is up to little Sidney, mostly, and $s(t)$ is equal to $D(t)$ on Thursdays and zero every other day.

Code that up and we have a system dynamics model.

In [2]:
# what would make this code really cool is if I could use the "with" syntax on a pandas.DataFrame

# that requires a context manager

from contextlib import contextmanager

@contextmanager
def columns_in(df):
col_names = df.columns
col_list = [df[col] for col in col_names]
try:
yield tuple(col_list)
finally:
for i, col in enumerate(col_names):
df[col] = col_list[i]

In [3]:
T = 120
state = pandas.DataFrame(index=range(-7,T), columns=['C', 'D', 'q', 'r', 's'])

state.ix[-7:0, ['C']] = 70
state.ix[-7:0, ['D', 'q', 'r', 's']] = 0

# simulate
with columns_in(state) as (C,D,q,r,s):
for t in range(T-1):
C[t+1] = C[t] + q[t] - r[t]
if C[t+1] < 0:
C[t+1] = 0

D[t+1] = D[t] + r[t] - s[t]

q[t+1] = s[t+1-7] if t % 7 == 6 else 0
# special case, first pickup, need to deliver something
if t == 6:
q[t+1] = 70.

r[t+1] = 8.
s[t+1] = D[t] if t % 7 == 6 else 0

t = state.filter(['C','D'])
t.columns = ['Clean', 'Dirty']
t.plot(linewidth=2, alpha=.9, figsize=(11.5, 4.5))
plot(arange(-7,T)[::7], state.ix[::7,'C'], ':o', mec='k', ms=7, color='k', linewidth=2, alpha=.7, label='Weekly minimum')

axis([-5,60,-10,120])
grid()
legend()
xlabel('Time (Days)')
ylabel('Diapers')

Out[3]:
<matplotlib.text.Text at 0xa3cdeac>

Now I can simulate my nightmare scenario, forgetting to take the dirty diapers out one week

In [4]:
T = 120
state = pandas.DataFrame(index=range(-7,T), columns=['C', 'D', 'q', 'r', 's'])

state.ix[-7:0, 'C'] = 70
state.ix[-7:0, ['D', 'q', 'r', 's']] = 0

# simulate
with columns_in(state) as (C,D,q,r,s):
for t in range(T-1):
C[t+1] = C[t] + q[t] - r[t]
if C[t+1] < 0:
C[t+1] = 0

D[t+1] = D[t] + r[t] - s[t]

q[t+1] = s[t+1-7] if t % 7 == 6 else 0
# special case, first pickup, need to deliver something
if t == 6:
q[t+1] = 70.

r[t+1] = 8.
s[t+1] = D[t] if t % 7 == 6 else 0

# special case, forget to take the diapers out
if t == 34:
s[t+1] = 0

t = state.filter(['C','D'])
t.columns = ['Clean', 'Dirty']
t.plot(linewidth=2, alpha=.9, figsize=(11.5, 4.5))
plot(arange(-7,T)[::7], state.ix[::7,'C'], ':o', mec='k', ms=7, color='k', linewidth=2, alpha=.7, label='Weekly minimum')

axis([-5,100,-10,130])
grid()
legend()
xlabel('Time (Days)')
ylabel('Diapers')

Out[4]:
<matplotlib.text.Text at 0xa55c28c>
In []: