Introduction

System dynamics models provide a means for exploring complex systems and effects that can arise from concepts such as feedback loops and time delays. This type of modeling has applications in a wide variety of fields including biology, economics, operations management, and social sciences. Industry standard tools for implementing this modeling process include Vensim and AnyLogic, but require a budget and lack the ability to construct programmatically from within Python. Existing Python-based libraries such as PySD provide the means to run models created in other tools but not to build them directly.

Bayesian inference is a statistical tool for modeling with uncertainty and updating probability distributions based on potentially limited amounts of data. PyMC is an established library in the Python ecosystem that provides algorithms for Bayesian inference, but it can be challenging to use for implementing complex system dynamics models. The goal of this project is to provide a Python-based means for building system dynamics models with a straightforward API, and support refinement of unknown or highly uncertain variables through PyMC without requiring the developer to write extensive PyMC specific code.

System Dynamics Implementation

We present Reno, a new open-source library with an API that centers around symbolically constructing equations that are used to define and reference stock, flow, and variable components, collectively constituting a system dynamics model. Conceptually similar to libraries like PyTensor and PyTorch, these equations create a compute graph that can be populated and evaluated to produce simulation data. Reno models, once defined, are called like a normal Python function to run a simulation, optionally passing in parameters to configure specific system variables. These model calls can efficiently run many simulations in parallel, allowing exploration of parameter space with parameter sweeps or input distributions, with results returned as XArray datasets.

This section will discuss an example from a system dynamics textbook and show the process of implementing it in Reno along with possible visualizations and analyses of the system once created.

Incorporating Bayesian Inference

By default, a Reno equation evaluates by running corresponding NumPy operations on the data passing through the compute graph. Given the similar API of PyTensor, the mathematics library underlying PyMC, everything within Reno compute graphs can also directly translate into a set of PyTensor/PyMC operations. A Reno model is thus converted into a PyMC model by compiling the component equations that evaluate for a single timestep, then wrapping with the necessary boilerplate to initialize the model and run the timestep function for a full time series simulation. Reno encapsulates this conversion process with a single function call, requiring no additional PyMC code from the model developer. Any observed data or measurements that are included in the function call are set within likelihood distributions and subsequently used in PyMC's MCMC sampling algorithms to approximate posterior distributions.

This section will expand on the previous example, showing how an uncertain input variable can be provided a prior probability distribution to indicate incomplete or imperfect knowledge. Further demonstration will show how the distribution tightens/converges around the ground truth value as additional observed data points are supplied to the PyMC model calls.

Links

Project repository: https://github.com/ornl/reno
Example of a previous SciPy talk: https://youtu.be/uyfIQEoZPOo