Systems Dynamics in Julia?

Spending a lot of time on the winter break reading Donella Meadows’ books (Thinking in Systems, Limits to Growth and sequels). Seems like a natural extension of the DiffEq.jl framework. Has anyone started anything along the lines of stock-and-flow modeling?

I’m familiar with some of the work that Chris R and others have done modeling COVID-19 SEIR models, which can be seen as a subset of stock-and-flow models.

I’m particularly interested in macros such as were developed for DiffEqBiological.jl or Catalyst.jl to define reaction networks, since similar descriptions are likely to be useful in Systems Programming. One of the major tools used in this field is Vensim, which has a nice graphical interface to setting up complicated reaction networks of this sort. Seems like one could replace the graphical networks with a macro-generated DSL that could be similarly agile and flexible and could leverage the DiffEq.jl solvers, Pluto notebooks, etc.

If anyone knows of existing efforts along these lines, please let me know.

Thanks, and happy new year!


Classic book. That reminds me that I need to go find my old copy.

I know Modelica has had a few systems dynamics implementations over the years. So ModelingToolkit in principal could do something similar since a lot of systems dynamics is naturally posed as DAEs.

1 Like

You probably already know that there’s a free online version of Limits to Growth.

Thanks for your comments.

There is a Systems Dynamics package for Modelica, developed by Francois Cellier (I think). Since ModelingToolkit shares some ideas with Modelica (but is more general), perhaps it is worthwhile to look at Cellier’s package. See, e.g., Microsoft Word - modelica2008_world3.doc (

Vensim does causal (signal-flow) modeling, like Simulink.
You may find what you want in zekeriyasari/Causal.jl: Causal.jl - A modeling and simulation framework adopting causal modeling approach. (

Looks excellent! Thanks for your reply.

Ahh yes, systems dynamics is something I’ve been looking at for a decade. Actually, fun fact, a long time ago I wrote pieces for a scripting language for a system modeling software (see the “Chris Rackauckas, June 2012” at the bottom :wink:) which is now known as Numerus. So I’ve built and played with these things quite a bit. One of the big choices one has to make in this kind of model is how continuous is interpreted between components. There are really two choices: either every continuous piece is solved together, or in isolation as a cosimulation.

Most simulation tools take the later choice, with Causal.jl being one of them. This can be very natural and allows any simulation to be composed together. For example, piecing a stochastic differential equation with a delay differential equation requires you specify two sets of equations, two solvers, and link them. There are mathematical and numerical problems with this of course though. For one, the composition of those two equations actually falls into the class of stochastic delay differential equations, and so you need specialized solvers to actually accurately capture the composition of the two. The delay differential equation piece in particular will be overly confident in its steps and have O(dt^{1/2}) accuracy while it things it has a lot more (usually O(dt^5)), which is mostly a theoretical concern (because “the equations aren’t exact/real”) until you start looking into things like “why won’t this stay positive?” with the answer being induced numerical error. The other problem associated with this approach is the handling of algebraic loops where, if you want to keep separate simulations intact, you need to jump through some hoops to get stability, or sometimes it might not be possible. Indeed, these explicit systems can become differential-algebraic equations through the connections, sometimes higher index DAEs, so things like Pantelides etc. are required to simulate them efficiently and accurately, and breaking of algebraic loops is a crude approximation to this.

This is not to harp on Causal.jl though. For such multi-clocked simulation environments, I would say that Causal.jl is probably the best one that I’ve ever seen. It’s just that the approach, while being by far the most common, has some major limitations for someone looking to achieve top level performance.

In contrast, the ModelingToolkit approach, based on the Modelica approach, considers the whole system as a single unit and simulates it together. Handling of algebraic loops becomes the symbolic tooling around differential-algebraic equations. A global symbolic view of the set of equations allows for symbolic simplifications, index reduction, tearing, etc. to all happen automatically. Everything generates into a single simulation code that can be orders of magnitude faster than having separate simulators within the components. Numerical accuracy is kept because you are forced to use a solver appropriate for the whole. There is a downside that supporting some more advanced functionality, like stochastic pieces + delays, is much more difficult, and sometimes the “faked” composed version is “fine enough” (and if it is fine enough for someone, it can be a lot faster than resorting to full SDE simulations which are much more expensive). However, I plan to get ModelingToolkit to cover those cases anyways, so we’ll get there.

Thus expect to see systems dynamics examples coming out of ModelingToolkit and related graphical tooling of JuliaSim.


Thanks for the reply, Chris.

Spent some time looking through Causal.jl today. Not exactly what I’m looking for. I can already model system dynamics perfectly well using DifferentialEquations.jl. But I’d like a compact form to describe the system.

I’d like something along the lines of the first choice you outline.
Catalyst.jl does almost everything I want. I haven’t looked through the Catalyst code (much), but the way that programs like Chemkin and Cantera work is that you have a reaction mechanism that is used to generate a Jacobian that is integrated, and I’m assuming that’s what you do in Catalyst.

It’s relatively easy to describe a SEIR model using what you already have in Catalyst:

rs = @reaction_network begin
  k1, S + I --> E + I
  k2, E --> I
  k3, I --> R
end k1 k2 k3

Maybe it’s just because I’ve stared at reaction networks for a long time, but there’s a lot of intuition that this model captures that might even be simpler than the graphical approaches.

I haven’t looked through something like World3 to see whether everything there can be described in elementary reactions. The parameters in those models get so convoluted I’m not convinced.

Catalyst.jl is pretty cool. Here’s the SEIR model from the QuantEcon paper using Catalyst.

1 Like

Unfortunately, I have not yet really found time to play around with Julia coming from Wolfram Mathematica und lately their Modelica platform, SystemModeler. Coming (originally) from Vensim DSS and System Dynamics in my professional work, I was interested in Cellier’s System Dynamics library—and found it unfortunately very cumbersome (compared with Vensim).

While Cellier’s library gives you the object-oriented advantages of Modelica it is completely stuck in a signal-flow-only framework, e.g., using blocks. Thus you cannot take advantage of the (acausal) physical connectors in Modelica where flows are automatically added (in the System dynamics library you would have to have a stock for 2 incoming flows and another one for 3 incoming flows …).

To make a long story short: I have written an open source Modelica library of my own called Business Simulation that uses acausal connectors for System Dynamics modeling. Take a look and maybe using OpenModelica or eventually Modia.jl will pave the road for further adaptation?

Here is a simple SIR model done with the library:

BSL Documentation

Business Simulation on GitHub


Dear @rpmuller, my colleagues and I are contributors to the PySD library, which we use to translate Vensim and Xmile models to Python and run them. However, what’s relevant to this conversation, is that we are planning to extend the functionality of PySD to also be able to translate the models to Julia.

During translation process we create an abstract representation of the model consisting of Python data classes, which we then translate to a string of Python code that we write to a file.

The advantage of having this abstract representation of the model, is that you can later use it to build the model in any other language of your choice (as long as you have written a model builder for that particular language).

That regards the translation part. Then we also need to build the code to import and solve the resulting model, which would be a separate library from PySD.

You can find more information on the paper we about to submit to the JOSS journal.

We hope to start working on designing the Julia builder inside the PySD library soon (Python code that writes Julia code), and also the solver (a standalone Julia package).

If anybody is interested to participate, even if it’s only a brainstorming, you will be more than welcome!


Very exciting package. Thanks for your work!