Getting Started with ModelingToolkit.jl

This is an introductory tutorial for ModelingToolkit (MTK). We will demonstrate the basics of the package by demonstrating how to define and simulate simple Ordinary Differential Equation (ODE) systems.

Installing ModelingToolkit

To install ModelingToolkit, use the Julia package manager. This can be done as follows:

using Pkg
Pkg.add("ModelingToolkit")

Copy-Pastable Simplified Example

A much deeper tutorial with forcing functions and sparse Jacobians is below. But if you want to just see some code and run, here's an example:

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

@mtkmodel FOL begin
    @parameters begin
        τ # parameters
    end
    @variables begin
        x(t) # dependent variables
    end
    @equations begin
        D(x) ~ (1 - x) / τ
    end
end

using DifferentialEquations: solve
@mtkbuild fol = FOL()
prob = ODEProblem(fol, [fol.x => 0.0], (0.0, 10.0), [fol.τ => 3.0])
sol = solve(prob)

using Plots
plot(sol)
Example block output

Now let's start digging into MTK!

Your very first ODE

Let us start with a minimal example. The system to be modelled is a first-order lag element:

\[\dot{x} = \frac{f(t) - x(t)}{\tau}\]

Here, $t$ is the independent variable (time), $x(t)$ is the (scalar) unknown variable, $f(t)$ is an external forcing function, and $\tau$ is a parameter. In MTK, this system can be modelled as follows. For simplicity, we first set the forcing function to a time-independent value $1$. And the independent variable $t$ is automatically added by $@mtkmodel$.

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

@mtkmodel FOL begin
    @parameters begin
        τ # parameters
    end
    @variables begin
        x(t) # dependent variables
    end
    @equations begin
        D(x) ~ (1 - x) / τ
    end
end

@mtkbuild fol = FOL()

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \frac{1 - x\left( t \right)}{\tau} \end{align} \]

Note that equations in MTK use the tilde character (~) as equality sign.

@mtkbuild creates an instance of FOL named as fol.

After construction of the ODE, you can solve it using DifferentialEquations.jl:

using DifferentialEquations
using Plots

prob = ODEProblem(fol, [fol.x => 0.0], (0.0, 10.0), [fol.τ => 3.0])
plot(solve(prob))
Example block output

The initial unknown and the parameter values are specified using a mapping from the actual symbolic elements to their values, represented as an array of Pairs, which are constructed using the => operator.

Algebraic relations and structural simplification

You could separate the calculation of the right-hand side, by introducing an intermediate variable RHS:

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

@mtkmodel FOL begin
    @parameters begin
        τ # parameters
    end
    @variables begin
        x(t) # dependent variables
        RHS(t)
    end
    @equations begin
        RHS ~ (1 - x) / τ
        D(x) ~ RHS
    end
end

@mtkbuild fol = FOL()

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \mathrm{RHS}\left( t \right) \end{align} \]

You can look at the equations by using the command equations:

equations(fol)

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \mathrm{RHS}\left( t \right) \end{align} \]

Notice that there is only one equation in this system, Differential(t)(x(t)) ~ RHS(t). The other equation was removed from the system and was transformed into an observed variable. Observed equations are variables which can be computed on-demand but are not necessary for the solution of the system, and thus MTK tracks it separately. One can check the observed equations via the observed function:

observed(fol)

\[ \begin{align} \mathrm{RHS}\left( t \right) =& \frac{1 - x\left( t \right)}{\tau} \end{align} \]

For more information on this process, see Observables and Variable Elimination.

MTK still knows how to calculate them out of the information available in a simulation result. The intermediate variable RHS therefore can be plotted along with the unknown variable. Note that this has to be requested explicitly like as follows:

prob = ODEProblem(fol,
    [fol.x => 0.0],
    (0.0, 10.0),
    [fol.τ => 3.0])
sol = solve(prob)
plot(sol, vars = [fol.x, fol.RHS])
Example block output

Named Indexing of Solutions

Note that the indexing of the solution similarly works via the names, and so to get the time series for x, one would do:

sol[fol.x]
14-element Vector{Float64}:
 0.0
 3.333277778395056e-5
 0.00036659945265974064
 0.0036931634343634343
 0.03635598665335361
 0.13698489268755668
 0.28615114066001374
 0.44421130189743585
 0.5989648996141732
 0.7323994349297617
 0.8366363895140377
 0.9093330812230745
 0.9548191028045532
 0.9643240555350271

or to get the second value in the time series for x:

sol[fol.x, 2]
3.333277778395056e-5

Similarly, the time series for RHS can be retrieved using the same indexing:

sol[fol.RHS]
14-element Vector{Float64}:
 0.3333333333333333
 0.33332222240740533
 0.3332111335157801
 0.3321022788552122
 0.3212146711155488
 0.28767170243748114
 0.23794961977999543
 0.18526289936752138
 0.1336783667952756
 0.0892001883567461
 0.05445453682865409
 0.03022230625897515
 0.015060299065148941
 0.011891981488324302

Specifying a time-variable forcing function

What if the forcing function (the “external input”) $f(t)$ is not constant? Obviously, one could use an explicit, symbolic function of time:

@mtkmodel FOL begin
    @parameters begin
        τ # parameters
    end
    @variables begin
        x(t) # dependent variables
        f(t)
    end
    @equations begin
        f ~ sin(t)
        D(x) ~ (f - x) / τ
    end
end

@named fol_variable_f = FOL()

\[ \begin{align} f\left( t \right) =& \sin\left( t \right) \\ \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \frac{ - x\left( t \right) + f\left( t \right)}{\tau} \end{align} \]

But often this function might not be available in an explicit form. Instead the function might be provided as time-series data. MTK handles this situation by allowing us to “register” arbitrary Julia functions, which are excluded from symbolic transformations, and thus used as-is. So, you could, for example, interpolate given the time-series using DataInterpolations.jl. Here, we illustrate this option by a simple lookup ("zero-order hold") of a vector of random values:

value_vector = randn(10)
f_fun(t) = t >= 10 ? value_vector[end] : value_vector[Int(floor(t)) + 1]
@register_symbolic f_fun(t)

@mtkmodel FOLExternalFunction begin
    @parameters begin
        τ # parameters
    end
    @variables begin
        x(t) # dependent variables
        f(t)
    end
    @structural_parameters begin
        h = 1
    end
    @equations begin
        f ~ f_fun(t)
        D(x) ~ (f - x) / τ
    end
end

@mtkbuild fol_external_f = FOLExternalFunction()
prob = ODEProblem(fol_external_f,
    [fol_external_f.x => 0.0],
    (0.0, 10.0),
    [fol_external_f.τ => 0.75])

sol = solve(prob)
plot(sol, vars = [fol_external_f.x, fol_external_f.f])
Example block output

Building component-based, hierarchical models

Working with simple one-equation systems is already fun, but composing more complex systems from simple ones is even more fun. Best practice for such a “modeling framework” could be to use factory functions for model components:

function fol_factory(separate = false; name)
    @parameters τ
    @variables x(t) f(t) RHS(t)

    eqs = separate ? [RHS ~ (f - x) / τ,
        D(x) ~ RHS] :
          D(x) ~ (f - x) / τ

    ODESystem(eqs, t; name)
end
fol_factory (generic function with 2 methods)

Such a factory can then be used to instantiate the same component multiple times, but allows for customization:

@named fol_1 = fol_factory()
@named fol_2 = fol_factory(true) # has observable RHS

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \mathrm{RHS}\left( t \right) \\ \mathrm{RHS}\left( t \right) =& \frac{ - x\left( t \right) + f\left( t \right)}{\tau} \end{align} \]

The @named macro rewrites fol_2 = fol_factory(true) into fol_2 = fol_factory(true,:fol_2). Now, these two components can be used as subsystems of a parent system, i.e. one level higher in the model hierarchy. The connections between the components again are just algebraic relations:

connections = [fol_1.f ~ 1.5,
    fol_2.f ~ fol_1.x]

connected = compose(ODESystem(connections, t, name = :connected), fol_1, fol_2)

\[ \begin{align} \mathrm{fol}_{1_{+}f}\left( t \right) =& 1.5 \\ \mathrm{fol}_{2_{+}f}\left( t \right) =& \mathrm{fol}_{1_{+}x}\left( t \right) \\ \frac{\mathrm{d} \mathrm{fol}_{1_{+}x}\left( t \right)}{\mathrm{d}t} =& \frac{\mathrm{fol}_{1_{+}f}\left( t \right) - \mathrm{fol}_{1_{+}x}\left( t \right)}{fol_{1_+\tau}} \\ \frac{\mathrm{d} \mathrm{fol}_{2_{+}x}\left( t \right)}{\mathrm{d}t} =& \mathrm{fol}_{2_{+}RHS}\left( t \right) \\ \mathrm{fol}_{2_{+}RHS}\left( t \right) =& \frac{\mathrm{fol}_{2_{+}f}\left( t \right) - \mathrm{fol}_{2_{+}x}\left( t \right)}{fol_{2_+\tau}} \end{align} \]

All equations, variables, and parameters are collected, but the structure of the hierarchical model is still preserved. This means you can still get information about fol_1 by addressing it by connected.fol_1, or its parameter by connected.fol_1.τ. Before simulation, we again eliminate the algebraic variables and connection equations from the system using structural simplification:

connected_simp = structural_simplify(connected)

\[ \begin{align} \frac{\mathrm{d} \mathrm{fol}_{1_{+}x}\left( t \right)}{\mathrm{d}t} =& \frac{\mathrm{fol}_{1_{+}f}\left( t \right) - \mathrm{fol}_{1_{+}x}\left( t \right)}{fol_{1_+\tau}} \\ \frac{\mathrm{d} \mathrm{fol}_{2_{+}x}\left( t \right)}{\mathrm{d}t} =& \mathrm{fol}_{2_{+}RHS}\left( t \right) \end{align} \]

full_equations(connected_simp)

\[ \begin{align} \frac{\mathrm{d} \mathrm{fol}_{1_{+}x}\left( t \right)}{\mathrm{d}t} =& \frac{1.5 - \mathrm{fol}_{1_{+}x}\left( t \right)}{fol_{1_+\tau}} \\ \frac{\mathrm{d} \mathrm{fol}_{2_{+}x}\left( t \right)}{\mathrm{d}t} =& \frac{\mathrm{fol}_{1_{+}x}\left( t \right) - \mathrm{fol}_{2_{+}x}\left( t \right)}{fol_{2_+\tau}} \end{align} \]

As expected, only the two equations with the derivatives of unknowns remain, as if you had manually eliminated as many variables as possible from the equations. Some observed variables are not expanded unless full_equations is used. As mentioned above, the hierarchical structure is preserved. So, the initial unknown and the parameter values can be specified accordingly when building the ODEProblem:

u0 = [fol_1.x => -0.5,
    fol_2.x => 1.0]

p = [fol_1.τ => 2.0,
    fol_2.τ => 4.0]

prob = ODEProblem(connected_simp, u0, (0.0, 10.0), p)
plot(solve(prob))
Example block output

More on this topic may be found in Composing Models and Building Reusable Components.

Initial Guess

It is often a good idea to specify reasonable values for the initial unknown and the parameters of a model component. Then, these do not have to be explicitly specified when constructing the ODEProblem.

@mtkmodel UnitstepFOLFactory begin
    @parameters begin
        τ = 1.0
    end
    @variables begin
        x(t) = 0.0
    end
    @equations begin
        D(x) ~ (1 - x) / τ
    end
end
ModelingToolkit.Model{typeof(Main.__UnitstepFOLFactory__), Dict{Symbol, Any}}(Main.__UnitstepFOLFactory__, Dict{Symbol, Any}(:variables => Dict{Symbol, Dict{Symbol, Any}}(:x => Dict(:default => 0.0, :type => Real)), :kwargs => Dict{Symbol, Dict}(:τ => Dict{Symbol, Any}(:value => 1.0, :type => Real), :x => Dict{Symbol, Any}(:value => 0.0, :type => Real)), :independent_variable => t, :parameters => Dict{Symbol, Dict{Symbol, Any}}(:τ => Dict(:default => 1.0, :type => Real)), :equations => Any["D(x) ~ (1 - x) / τ"]), false)

While defining the model UnitstepFOLFactory, an initial guess of 0.0 is assigned to x(t) and 1.0 to τ. Additionally, these initial guesses can be modified while creating instances of UnitstepFOLFactory by passing arguments.

@mtkbuild fol = UnitstepFOLFactory(; x = 0.1)
sol = ODEProblem(fol, [], (0.0, 5.0), []) |> solve
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 12-element Vector{Float64}:
 0.0
 0.0645677678571118
 0.23112586321494216
 0.47279997996782797
 0.7782099882421787
 1.168340018958489
 1.6418244200442873
 2.2128069066549645
 2.888015454160246
 3.6845362502170076
 4.6215123142148355
 5.0
u: 12-element Vector{Vector{Float64}}:
 [0.1]
 [0.15627467655510244]
 [0.285724385884508]
 [0.4390707362116239]
 [0.5866953750363533]
 [0.7202054759739022]
 [0.8257356363164883]
 [0.9015449468569201]
 [0.9498798747479449]
 [0.97739869398512]
 [0.9911413585329305]
 [0.9939327453112139]

In non-DSL definitions, one can pass defaults dictionary to set the initial guess of the symbolic variables.

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

function UnitstepFOLFactory(; name)
    @parameters τ
    @variables x(t)
    ODESystem(D(x) ~ (1 - x) / τ; name, defaults = Dict(x => 0.0, τ => 1.0))
end
UnitstepFOLFactory (generic function with 1 method)

Note that the defaults can be functions of the other variables, which is then resolved at the time of the problem construction. Of course, the factory function could accept additional arguments to optionally specify the initial unknown or parameter values, etc.

Symbolic and sparse derivatives

One advantage of a symbolic toolkit is that derivatives can be calculated explicitly, and that the incidence matrix of partial derivatives (the “sparsity pattern”) can also be explicitly derived. These two facts lead to a substantial speedup of all model calculations, e.g. when simulating a model over time using an ODE solver.

By default, analytical derivatives and sparse matrices, e.g. for the Jacobian, the matrix of first partial derivatives, are not used. Let's benchmark this (prob still is the problem using the connected_simp system above):

using BenchmarkTools
@btime solve(prob, Rodas4());
  46.006 μs (158 allocations: 18.11 KiB)

Now have MTK provide sparse, analytical derivatives to the solver. This has to be specified during the construction of the ODEProblem:

prob_an = ODEProblem(connected_simp, u0, (0.0, 10.0), p; jac = true)
@btime solve($prob_an, Rodas4());
  35.857 μs (119 allocations: 14.75 KiB)
prob_an = ODEProblem(connected_simp, u0, (0.0, 10.0), p; jac = true, sparse = true)
@btime solve($prob_an, Rodas4());
  121.987 μs (1147 allocations: 70.62 KiB)

The speedup is significant. For this small dense model (3 of 4 entries are populated), using sparse matrices is counterproductive in terms of required memory allocations. For large, hierarchically built models, which tend to be sparse, speedup and the reduction of memory allocation can be expected to be substantial. In addition, these problem builders allow for automatic parallelism using the structural information. For more information, see the ODESystem page.

Notes and pointers how to go on

Here are some notes that may be helpful during your initial steps with MTK:

  • The @mtkmodel macro is for high-level usage of MTK. However, in many cases you may need to programmatically generate ODESystems. If that's the case, check out the Programmatically Generating and Scripting ODESystems Tutorial.
  • Vector-valued parameters and variables are possible. A cleaner, more consistent treatment of these is still a work in progress, however. Once finished, this introductory tutorial will also cover this feature.

Where to go next?