Building a System Dynamics Model in Python with PySD

Last updated: 31 October 2026. Code tested with PySD 3.14.3 on Python 3.13.

System dynamics in Python lets you keep the strengths of a stock-and-flow model, a clear structure and explicit feedback, while using the Python ecosystem for data handling, batch runs and plotting. PySD is the open-source library most often used for this. It reads models built in Vensim or in the XMILE format that Stella and others export, translates them to Python, and runs them. This tutorial builds a small adoption model, runs it, checks the result against a known closed-form answer, and then uses parameter sweeps to test which assumptions matter.

What PySD does and when to use it

PySD (the name stands for Python System Dynamics) translates a model file into Python code, then lets you run and modify it through an object. According to its documentation it supports models from Vensim and from XMILE-compatible tools. The practical benefits:

• Batch experiments. Running hundreds of scenarios in a loop rather than clicking through a GUI.
• Data integration. Reading inputs from files or databases and merging results with other analyses.
• Reproducibility. A script and a model file are easy to version and re-run.
• Calibration and uncertainty analysis. You can wrap a model in optimisation or Monte Carlo code.

Use it when you already have a structural model and need more experiments or data plumbing than the modelling tool provides. If you are still deciding on structure, stay in a tool where you can draw and inspect stocks, flows and loops; see stocks and flows in system dynamics for the building blocks.

Install

pip install pysd

Use a virtual environment so the pandas and xarray versions that PySD depends on do not clash with other projects.

The model: product adoption with advertising and word of mouth

The structure is a classic. Potential adopters become adopters through two channels: outside influence such as advertising, and word of mouth from existing adopters meeting non-adopters. The second channel is a reinforcing loop, the first is a steady external push. This is the Bass diffusion model, rewritten as stocks and flows.

Figure 1. Left: the model structure. Right: the adopters stock simulated in PySD over 36 months, for the baseline and for a run with the word-of-mouth conversion probability halved.

Here is the model in Vensim text format. Save it as adoption.mdl (the equations are the substance; the control block sets a 36-month run with a time step of 0.0625 month).

Potential Adopters= INTEG(-Adoption Rate, Total Population-Initial Adopters)
    ~ people ~ |
Adopters= INTEG(Adoption Rate, Initial Adopters)
    ~ people ~ |
Adoption Rate= Advertising Adoption + Word of Mouth Adoption
    ~ people/Month ~ |
Advertising Adoption= Potential Adopters*Advertising Effectiveness
    ~ people/Month ~ |
Word of Mouth Adoption= Contact Rate*Adoption Fraction*Adopters*Potential Adopters/Total Population
    ~ people/Month ~ |
Advertising Effectiveness= 0.03 ~ 1/Month ~ |
Contact Rate= 20 ~ 1/Month ~ |
Adoption Fraction= 0.019 ~ Dmnl ~ |
Initial Adopters= 10 ~ people ~ |
Total Population= 10000 ~ people ~ |
FINAL TIME= 36 ~ Month ~ |
INITIAL TIME= 0 ~ Month ~ |
SAVEPER= TIME STEP ~ Month ~ |
TIME STEP= 0.0625 ~ Month ~ |

Two details of the equations are worth noting. The word-of-mouth term multiplies adopters by the fraction of the population still available (Potential Adopters / Total Population), which is the balancing part that makes the curve saturate. And the units check out: contacts per month times a dimensionless conversion probability times people gives people per month.

Load and run

import pysd

model = pysd.read_vensim("adoption.mdl")
results = model.run(return_columns=["Adopters", "Adoption Rate"])

print(results.head(3))
print(results["Adoption Rate"].idxmax(), round(results["Adoption Rate"].max()))

The first call translates the file into a Python model. The second runs the simulation using the model’s own control settings and returns a pandas DataFrame indexed by time. In my run the first rows were:

        Adopters  Adoption Rate
time
0.0000      10.0          303.5
0.0625      29.0          310.1
0.1250      48.4          316.8

and the adoption rate peaked at about 1,106 people per month at month 6.19. The model.doc attribute returns a table of every variable with its units, type (constant, auxiliary or stateful) and comment, which is the fastest way to check what PySD understood from the file.

Validate against a known answer

Before trusting any run, check it against something you can derive independently. For the Bass model the peak adoption rate occurs at time ln(q/p)/(p+q), with p the advertising coefficient (0.03) and q the word-of-mouth coefficient, here contact rate times conversion: 20 × 0.019 = 0.38. That gives ln(0.38/0.03)/0.41, about 6.19 months, and a peak rate of (p+q)²/(4q) times the population, about 1,106 per month. The simulation agrees. This kind of check is a minimal form of the practice described in model validation, and it catches errors such as wrong units or a mistyped equation that a plausible-looking curve would hide.

For a model with no closed form, use other tests: extreme conditions (set advertising to zero and adopters should still grow only if the initial adopters are above zero), mass conservation (potential adopters plus adopters should always equal the total population), and replicating a published result.

Change parameters and compare scenarios

for f in (0.0095, 0.019, 0.038):
    m = pysd.read_vensim("adoption.mdl")
    r = m.run(params={"Adoption Fraction": f}, return_columns=["Adopters"])
    month = r.index[(r["Adopters"] >= 5000).argmax()]
    print(f, round(float(month), 2))

The params argument overrides a constant for that run. Creating a fresh model object for each run keeps runs independent. The output:

Adoption Fraction	Months until half the market has adopted
0.0095 (half)	9.69
0.019 (baseline)	6.56
0.038 (double)	4.25

Halving the conversion probability delays the half-way point by about three months, and doubling it brings it forward by about two and a third. The response is not symmetric, which is typical of reinforcing loops, and a sweep like this shows it at once. For a more formal treatment see sensitivity analysis in systems models.

Plot the behaviour over time

Because the result is a DataFrame, plotting is one line with pandas and Matplotlib:

results["Adopters"].plot(ylabel="people")

The shape to expect is the S-curve in Figure 1. Reading such charts is a skill in itself; the guide to behaviour-over-time graphs covers what to look for.

The same model in plain Python

It is worth seeing how little a model of this size needs. The following loop uses Euler integration, as the PySD run above does, and needs no libraries:

N, A, P, dt, t = 10000, 10.0, 9990.0, 0.0625, 0.0
peak = (0, 0)
while t <= 36:
    rate = P*0.03 + 20*0.019*A*P/N
    if rate > peak[1]:
        peak = (t, rate)
    A += rate*dt
    P -= rate*dt
    t += dt
print(round(peak[0], 2), round(peak[1]))

It prints the same peak, 6.19 months and 1,106 people per month. This is a useful sanity check on PySD results, and a reminder that the library’s value is for larger models with many variables, units, lookups and delays, where writing and checking the equations by hand becomes error-prone.

Practical tips and pitfalls

• Check the translation. Use model.doc and test a few values. Not every Vensim function is supported by every PySD version, and the documentation lists the supported set.
• Mind the time step. Integration error depends on it. Halve the time step and confirm the results hardly change.
• One model object per independent run. Models keep state, so reuse across scenarios can leak state unless you reset it deliberately.
• Return only the columns you need. With large models, return_columns keeps memory and run time down.
• Keep the model file as the source of truth. Do not patch equations in Python and forget to update the diagram.
• Link results back to the question. Behaviour-over-time and sensitivity plots exist to answer a policy question, not to demonstrate the tool.

Make the model reproducible

A model that runs once is not yet a reliable analysis. Keep the source model, imported Python representation, input data and package versions together so a colleague can reproduce the run. Record the simulation start and end, time step, units and any initialization choices. These details matter because a mismatch in time units or an unnoticed change to a lookup table can create a plausible graph with the wrong meaning.

Build small checks around the model. Test that stocks remain nonnegative where the real quantity cannot fall below zero, that flows have the right units, and that a deliberately extreme parameter produces an understandable response. Compare one run with a known result, as this tutorial does, then test behavior across a range of inputs. For stochastic models, set and record the random seed while also checking results across multiple seeds.

Version control makes this workflow safer. Commit the original model and the scripts used to run it; save parameter sets as data rather than editing code by hand between scenarios. A scenario table can then show which assumption changed and what outcome moved. This supports review: a decision maker can distinguish a structural change to the model from a new value for an uncertain parameter.

Frequently asked questions

Does PySD work with Stella models?

PySD reads XMILE files, the interchange format that Stella and other tools can export, through pysd.read_xmile(). Check the documentation for which features your version supports.

Can I build a model from scratch in Python with PySD?

PySD is designed to translate and run existing models rather than as a drawing environment. For building, use a modelling tool, or write the equations directly in Python if the model is small, as in the plain Python version above.

How do I run many scenarios?

Loop over parameter values and pass them through params= on each run, as shown in the sweep. For larger designs, combine this with a sampling library and run in parallel processes.

How do I know the model is right?

You cannot prove it, but you can test it: compare with a closed-form case, check conservation, test extreme conditions, and compare behaviour with real data. The reference mode (the observed pattern the model should reproduce) is the starting point for the last of these.

Conclusion

PySD is a bridge between the structural clarity of system dynamics and the scripting power of Python. The workflow is short: build and inspect the structure in a modelling tool, load it with read_vensim, run and validate against any analytical case you can find, then use parameter sweeps to see which assumptions drive the result. The checks matter more than the tooling: a model that has not been validated is an illustration, not evidence.

Related terms: PySD, system dynamics, stock and flow. For emergent patterns you may want to simulate rather than solve, see weak versus strong emergence.

Sources and further reading

• PySD documentation, including model loading and running.
• Bass, F. M., “A new product growth model for consumer durables”, Management Science, 1969.
• Sterman, J. D., Business Dynamics: Systems Thinking and Modeling for a Complex World, McGraw-Hill, 2000.