Systems Modeling from scratch
I recently finished Donella H. Meadow’s “Thinking in Systems,” a primer on systems modeling. The book shows how phenomena like ecosystems, economies, business inventory, and room temperature can be modeled as a system where stocks of resources change at rates based on flows that depend on information. Understanding the structure will reveal reinforcing and balancing feedback loops and nonlinearities, and can provide strategies for changing systems.
One part of the book is about interesting and counterintuitive behavior that appears when running simulations for systems with certain structures. I wanted to try exploring some of these. The systems Meadows discuss are a graph with computations. As that lends itself to code, I wrote a 250-line simulator.
This post shows using the script to implement three example systems from the book. The systems and equations are transcribed from the book’s appendix, and more complex equations are made up to match the graphs.
Defining a system
Before diving into the example systems, this section describes how my script’s
Components map to those in System Dynamics. It then shows how a
Component is defined and running the simulation.
I tried to code
system_simulation.py to mostly match how Meadows describes systems.
Two components of the system are stocks (
StockComponent), which remember how much of something there is, and flows (
FlowComponent), which measure the rate of change of the stocks.
An example of a system with stocks and flows is a bathtub: The stock is the amount of water in a bathtub. The inflow is the rate at which the faucet pours water into the bathtub. The outflow is the rate at which the water leaves the bathtub through the drain. The total amount of water in the bathtub is
In these systems, stocks cannot change except through flows. The flow can change depending on a
StockComponent or a third component, the
InfoComponent, which represents how information travels between stocks, flows, and other intermediate information components. (In the appendix, Meadows calls
Defining a system
room_temperature depends on
heat_from_outside, and has an initial value of 10 degrees Celsius, I can define a stock with code like this:
StockComponent( 'room_temperature', initial_value=10, inflow='heat_from_furnace', outflow='heat_to_outside', min_value=-273.15 )
Then I’d need to define the inflow and outflow. The
FlowComponents can define an equation for the flow, which takes the current timestep
t and an argument for each of the
parents in order.
FlowComponent( 'heat_from_furnace', lambda t, thermostat_discrepancy: min(thermostat_discrepancy, 5), initial_value=0, parents=[Parent('thermostat_discrepancy')] ) FlowComponent( 'heat_to_outside', lambda t, outdoor_discrepancy: outdoor_discrepancy * .1, initial_value=0, parents=[Parent('outdoor_discrepancy')] )
The rest of the connections are called
InfoComponents, which look the same way as
- Normally the equations are passed the most recent value of the parent. To encode things like delays,
Parent(..., with_history=True)will pass the equation the entire history of that component.
- Loops are one of the most interesting parts of systems. But computations need to be done sequentially. To encode a feedback loop,
Parent(..., prev_timestep=True), will tell the equation to read the value computed at the previous timestep.
StockComponentby default use the in- and outflows computed in the previous timestep.
After defining the network, I run the simulations for some number of steps.
samples_per_hour = 2 hours = 100 t_num = hours * samples_per_hour dt = 1/samples_per_hour simulation = s.simulate(t_num, dt)
The above code would simulate 8 hours at a granularity of 30 minutes. It does so by computing values at 16 timesteps using a step size of 0.5 hours.
The simulation runs by computing each node for each timestep.
simulation contains the values of each node over the timesteps.
A few models
What follows are three models from the book “Thinking in Systems.” If you’d like to know more about the structure, equations, and resulting behavior, the book describes each system in detail in Chapter 2.
System 1: Thermostat
The Thermostat system represents a room on a cold day. A furnace warms the room and heat leaks from the room.
The thermostat system has two balancing (negative) feedback loops.
The results of the simulation are shown below. The room gradually warms up to just below the thermostat setting.
System 2: Inventory with delays
The next system models the inventory of a business after an increase in customer demand. This system has a similar structure to the room temperature system, except that delays are introduced.
- The delivery delay means that it takes a few days for orders to become deliveries.
- Because sales may have random noise, it’s important to soften reactions to them. This perception delay means that the perceived sales (which orders depend on) is the average of actual sales over a few timesteps.
- The response delay controls how quickly to try to make up the discrepancy between desired and actual inventory.
The inventory system also has one stock (inventory) and two balancing feedback loops.
Simulation with no delays
In the case of no delays, the inventory quickly adjusts in response to the change in customer demand.
The red line shows when the customer demand changed.
Simulation with delays
When delays are introduced, wild oscillations happen in response to a small change in demand.
Simulation with different response delays
The first reaction may be that “delay is bad, reduce all the delay”. The delivery delay is outside of our control, and the perception delay is a good idea to account for noise. But counterintuitively, reducing the response delay results in larger oscillations! Instead, increasing the response delay will lead to a smoother reaction.
System 3: Renewable Resource
The renewable resource system is a more complicated system from the book.
There are two stocks:
- the renewable resource of fish (stock), which regenerate (inflow) and are harvested (outflow)
- the capital used to catch more fish (stock), which grows with investment (inflow), and depreciates (outflow)
In addition, there are many other connections which introduce feedback loops. For example, there is a reinforcing feedback loop of capital, which creates more investment which increases capital and so on.
Aside, some helper functions
The appendix of the book plots a few magical functions. I tried to approximate them below.
regeneration_rate_given_resource, that gives some optimal regeneration rate location. If there are too few fish, they can’t find each other. If there are too many fish, they fight over food.
The third function
yield_per_unit_capital_given_resource, turns out to be the most important this example.
The example given in the book is how much money it takes to get the same number of fish. Moving the location of the curve to the right might represent boats that can go into deeper water and catch more fish with the same amount of capital.
The book shows how different curves result in different scenerios. I admit I don’t know the theory behind it, so I just call the parameter
some_measure_of_efficiency, which shifts the curve to the left or right.
If they are too good at catching fish, the number of fish drops to nearly 0.
Simulation: Equal out
If they are a little worse at catching fish, the number of fish stabilizes.
Aside: Fun plots
If instead of plotting resource and capital over time, I plot resource and capital, it shows a cool spiral as the values approach an attractor.
Simulation: Oscillations forever
If they are even worse at catching fish, the number of fish oscillates.
Aside: Fun plots
This post gave a quick tour of interesting behavior from some systems based on the example models from “Thinking in Systems” using a quick implementation
Some additional places of the web that might be interesting:
- Thinking in Systems
- Wikipedia on System dynamics
- For pro software for modeling system dynamics, there’s MATLAB’s Simulink