GitHub

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 it, here's an example:

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

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

using OrdinaryDiffEq
@mtkcompile fol = FOL()
prob = ODEProblem(fol, [], (0.0, 10.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
        τ = 3.0 # parameters and their values
    end
    @variables begin
        x(t) = 0.0 # dependent variables and their initial conditions
    end
    @equations begin
        D(x) ~ (1 - x) / τ
    end
end

@mtkcompile 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.

@mtkcompile creates an instance of FOL named as fol.

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

using OrdinaryDiffEq
using Plots

prob = ODEProblem(fol, [], (0.0, 10.0))
plot(solve(prob))
Example block output

The parameter values are determined using the right hand side of the expressions in the @parameters block, and similarly initial conditions are determined using the right hand side of the expressions in the @variables block.

Using different values for parameters and initial conditions

If you want to simulate the same model, but with different values for the parameters and initial conditions than the default values, you likely do not want to write an entirely new @mtkmodel. ModelingToolkit supports overwriting the default values:

@mtkcompile fol_different_values = FOL(; τ = 1 / 3, x = 0.5)
prob = ODEProblem(fol_different_values, [], (0.0, 10.0))
plot(solve(prob))
Example block output

Alternatively, this overwriting could also have occurred at the ODEProblem level.

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

Here, the second argument of ODEProblem is an array of Pairs. The left hand side of each Pair is the parameter you want to overwrite, and the right hand side is the value to overwrite it with. Similarly, the initial conditions are overwritten in the fourth argument. One important difference with the previous method is that the parameter has to be referred to as fol.τ instead of just τ.

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
        τ = 3.0 # parameters and their values
    end
    @variables begin
        x(t) = 0.0 # dependent variables and their initial conditions
        RHS(t)
    end
    @equations begin
        RHS ~ (1 - x) / τ
        D(x) ~ RHS
    end
end

@mtkcompile fol = FOL()

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

If you copy this block of code to your REPL, you will not see the above LaTeX equations. Instead, you can look at the equations by using the equations function:

equations(fol)

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} &= \mathtt{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 that can be computed on-demand but are not necessary for the solution of the system, and thus MTK tracks them separately. For this reason, we also did not need to specify an initial condition for RHS. You can check the observed equations via the observed function:

observed(fol)

\[ \begin{align} \mathtt{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:

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

Named Indexing of Solutions

Note that the indexing of the solution also works via the symbol, and so to get the time series for x, you would do:

sol[fol.x]
14-element Vector{Float64}:
 0.0
 3.333277778395056e-5
 0.00036659945265974064
 0.0036931634343634343
 0.03635598665335361
 0.1369848912590037
 0.2861511310403561
 0.4442112892126739
 0.5989648889011088
 0.7323994289466121
 0.8366363855665494
 0.9093330787024341
 0.9548191013493363
 0.9643240555351114

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 symbolic indexing:

sol[fol.RHS]
14-element Vector{Float64}:
 0.3333333333333333
 0.33332222240740533
 0.3332111335157801
 0.3321022788552122
 0.3212146711155488
 0.28767170291366545
 0.23794962298654795
 0.18526290359577535
 0.13367837036629707
 0.08920019035112929
 0.05445453814448353
 0.030222307099188644
 0.01506029955022122
 0.011891981488296214

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
        τ = 3.0 # parameters and their values
    end
    @variables begin
        x(t) = 0.0 # dependent variables and their initial conditions
        f(t)
    end
    @equations begin
        f ~ sin(t)
        D(x) ~ (f - x) / τ
    end
end

@mtkcompile fol_variable_f = FOL()

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

However, 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. For example, you could interpolate given the time-series using DataInterpolations.jl. Here, we illustrate this option with 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
        τ = 0.75 # parameters and their values
    end
    @variables begin
        x(t) = 0.0 # dependent variables and their initial conditions
        f(t)
    end
    @equations begin
        f ~ f_fun(t)
        D(x) ~ (f - x) / τ
    end
end

@mtkcompile fol_external_f = FOLExternalFunction()

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

prob = ODEProblem(fol_external_f, [], (0.0, 10.0))
sol = solve(prob)
plot(sol, idxs = [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. The best practice for such a “modeling framework” is to use the @components block in the @mtkmodel macro:

@mtkmodel FOLUnconnectedFunction begin
    @parameters begin
        τ # parameters
    end
    @variables begin
        x(t) # dependent variables
        f(t)
        RHS(t)
    end
    @equations begin
        RHS ~ f
        D(x) ~ (RHS - x) / τ
    end
end
@mtkmodel FOLConnected begin
    @components begin
        fol_1 = FOLUnconnectedFunction(; τ = 2.0, x = -0.5)
        fol_2 = FOLUnconnectedFunction(; τ = 4.0, x = 1.0)
    end
    @equations begin
        fol_1.f ~ 1.5
        fol_2.f ~ fol_1.x
    end
end
@mtkcompile connected = FOLConnected()

\[ \begin{align} \frac{\mathrm{d} \mathtt{fol\_2.x}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{fol\_2.RHS}\left( t \right) - \mathtt{fol\_2.x}\left( t \right)}{\mathtt{fol\_2.\tau}} \\ \frac{\mathrm{d} \mathtt{fol\_1.x}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{fol\_1.RHS}\left( t \right) - \mathtt{fol\_1.x}\left( t \right)}{\mathtt{fol\_1.\tau}} \end{align} \]

Here the total model consists of two of the same submodels (components), but with a different input function, parameter values and initial conditions. The first model has a constant input, and the second model uses the state x of the first system as an input. To avoid having to type the same differential equation multiple times, we define the submodel in a separate @mtkmodel. We then reuse this submodel twice in the total model @components block. The inputs of two submodels then still have to be specified in the @equations block.

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.τ.

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:

prob = ODEProblem(connected, [], (0.0, 10.0))
plot(solve(prob))
Example block output

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

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 system above):

using BenchmarkTools
@btime solve(prob, Rodas4());
  75.360 μs (254 allocations: 23.00 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, [], (0.0, 10.0); jac = true)
@btime solve(prob_an, Rodas4());
  65.943 μs (243 allocations: 20.09 KiB)
prob_sparse = ODEProblem(connected, [], (0.0, 10.0); jac = true, sparse = true)
@btime solve(prob_sparse, Rodas4());
  192.880 μs (1392 allocations: 82.49 KiB)

The speedup using the analytical Jacobian is significant. For this small dense model (3 of 4 entries 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 also be expected to be substantial. In addition, these problem builders allow for automatic parallelism by exploiting the structural information. For more information, see the System 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 Systems. If that's the case, check out the Programmatically Generating and Scripting Systems 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?