ModelingToolkit.jl

ModelingToolkit.jl is a modeling framework for high-performance symbolic-numeric computation in scientific computing and scientific machine learning. It allows for users to give a high-level description of a model for symbolic preprocessing to analyze and enhance the model. ModelingToolkit can automatically generate fast functions for model components like Jacobians and Hessians, along with automatically sparsifying and parallelizing the computations. Automatic transformations, such as index reduction, can be applied to the model to make it easier for numerical solvers to handle.

For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation which contains the unreleased features.

Standard Library

For a standard library of ModelingToolkit components and blocks, check out the ModelingToolkitStandardLibrary

High-Level Examples

First, let's define a second order riff on the Lorenz equations, symbolically lower it to a first order system, symbolically generate the Jacobian function for the numerical integrator, and solve it.

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

# Defines a ModelingToolkit `System` model.
@parameters σ ρ β
@variables x(t) y(t) z(t)
eqs = [
    D(D(x)) ~ σ * (y - x),
    D(y) ~ x * (ρ - z) - y,
    D(z) ~ x * y - β * z
]
@mtkcompile sys = System(eqs, t)

# Simulate the model for a specific condition (initial condition and parameter values).
using OrdinaryDiffEqDefault
sim_cond = [
    D(x) => 2.0,
    x => 1.0,
    y => 0.0,
    z => 0.0,
    σ => 28.0,
    ρ => 10.0,
    β => 8 / 3
]
tend = 100.0
prob = ODEProblem(sys, sim_cond, tend; jac = true)
sol = solve(prob)

# Plot the solution in phase-space.
using Plots
plot(sol, idxs = (x, y))

This will have automatically generated fast Jacobian functions, making it more optimized than directly building a function. In addition, we can then use ModelingToolkit to compose multiple ODE subsystems. Now, let's define two interacting Lorenz equations and simulate the resulting Differential-Algebraic Equation (DAE):

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

# Defines two lorenz system models.
eqs = [
    D(x) ~ σ * (y - x),
    D(y) ~ x * (ρ - z) - y,
    D(z) ~ x * y - β * z
]
@named lorenz1 = System(eqs, t)
@named lorenz2 = System(eqs, t)

# Connect the two models, creating a single model.
@variables a(t)
@parameters γ
connections = [0 ~ lorenz1.x + lorenz2.y + a * γ]
@mtkcompile connected_lorenz = System(connections, t; systems = [lorenz1, lorenz2])

# Simulate the model for a specific condition (initial condition and parameter values).
using OrdinaryDiffEqDefault
sim_cond = [
    lorenz1.x => 1.0,
    lorenz1.y => 0.0,
    lorenz1.z => 0.0,
    lorenz2.x => 0.0,
    lorenz2.z => 0.0,
    a => 2.0,
    lorenz1.σ => 10.0,
    lorenz1.ρ => 28.0,
    lorenz1.β => 8 / 3,
    lorenz2.σ => 10.0,
    lorenz2.ρ => 28.0,
    lorenz2.β => 8 / 3,
    γ => 2.0
]
tend = 100.0
prob = ODEProblem(connected_lorenz, sim_cond, tend)
sol = solve(prob)

# Plot the solution in phase-space.
using Plots
plot(sol, idxs = (a, lorenz1.x, lorenz2.z))

Citation

If you use ModelingToolkit.jl in your research, please cite this paper:

@misc{ma2021modelingtoolkit,
      title={ModelingToolkit: A Composable Graph Transformation System For Equation-Based Modeling},
      author={Yingbo Ma and Shashi Gowda and Ranjan Anantharaman and Chris Laughman and Viral Shah and Chris Rackauckas},
      year={2021},
      eprint={2103.05244},
      archivePrefix={arXiv},
      primaryClass={cs.MS}
}