ModelingToolkit.jl is a Julia package which provides a symbolic interface for specifying numerical problems such as systems of differential equations, optimization problems, optimal control problems and more. It is able to symbolically simplify such systems and generate efficient code to target numerical solvers.

At JuliaCon 2025, Aayush Sabharwal presented a talk about the growth of ModelingToolkit over the past couple of years, highlighting new features, improved workflows and documentation and plans for the future.

Feature Retrospective

MTKParameters

One of the most significant changes in ModelingToolkit version 9 was the introduction of the MTKParameters object. This is used as the parameter object by default for all MTK-generated code. It enables storing parameters of different types in a type-stable manner, where the previous approach would resort to a non-concrete array. For example, MTK can accept functions as parameters:

@variables x(t)
@parameters f(::Real)::Real q
@mtkcompile sys = System([D(x) ~ f(t) + q * x], t)
prob = ODEProblem(sys,
    [x => 1.0, f => sin, q => -1.0], (0.0, 5.0))
sol = solve(prob, Tsit5())
plot(sol)

Here, f is a parameter representing a function that takes a real number and outputs a real number. It is used in the equations as f(t). In the ODEProblem constructor, it is assigned the sin function. This generates the following result:

The value of f can be swapped out for any other function with the same interface, in the same way any other parameter value is changed. For example, it can be changed to represent the cos function:

prob.ps[f] = cos
sol = solve(prob, Tsit5())
plot(sol)

This now gives the following result:

MTKParameters also leverages the SciMLStructures.jl interface to enable efficient parameter optimization. Below, we create a simple system of ODEs with four parameters. The first three are tunable by default (meaning we want to optimize them) and the fourth explicitly opts out of this category:

@parameters α β γ δ [tunable = false]
@variables x(t) y(t)
eqs = [D(x) ~ (α - β * y) * x
       D(y) ~ (δ * x - γ) * y]
@mtkcompile odesys = System(eqs, t)

We can set up our loss function to optimize for the parameters. This uses the Optimization.jl interface, where x represents the current state of the three ODE parameters and p is the parameter object of the optimization.

function loss(x, p)
    odeprob = p[1] # ODEProblem stored as parameters to avoid using global variables
    ps = parameter_values(odeprob) # obtain the parameter object from the problem
    # create a copy of the parameter object with the buffer
    ps = SciMLStructures.replace(SciMLStructures.Tunable(), ps, x)
    # remake the problem, passing in our new parameter object
    newprob = remake(odeprob; p = ps)
    timesteps = p[2]
    sol = solve(newprob, AutoTsit5(Rosenbrock23()); saveat = timesteps)
    truth = p[3]
    data = Array(sol)
    return sum((truth .- data) .^ 2) / length(truth)
end

Note how the loss function takes the ODEProblem from p, obtains the parameter object from it, and calls SciMLStructures.replace(SciMLStructures.Tunable(), ps, x). This effectively means “replace the tunable part of ps with x”. The ODEProblem is then remade with the new parameter object, solved and the loss calculated against experimental data.

# experimental data to optimize against
timesteps = # ...
data = # ...
optfn = OptimizationFunction(loss, Optimization.AutoForwardDiff())
# parameter object is a tuple, to store differently typed objects together
optprob = OptimizationProblem(
    optfn, rand(3), (odeprob, timesteps, data),
    lb = 0.1zeros(3), ub = 3ones(3))
sol = solve(optprob, BFGS())
ps = parameter_values(odeprob)
ps = replace(Tunable(), ps, sol.u)
newprob = remake(odeprob; p = ps)
newprob.ps[[α, β, γ

Assuming the presence of experimental data, we can set up and solve the OptimizationProblem. The values from the resultant solution can be replaced into the ODEProblem in much the same way as in the loss function, and the new parameter values can be accessed symbolically from there.

Initialization

Initialization allows ModelingToolkit.jl to solve for the initial conditions of your system independent of which variables are selected as states. Take, for example, the classic pendulum:

@parameters g
@variables x(t)
@variables y(t) [state_priority = 10, guess = 1.0]
@variables λ(t) [guess = 1.0]
eqs = [D(D(x)) ~ λ * x
       D(D(y)) ~ λ * y - g
       x^2 + y^2 ~ 1]
@mtkcompile pend = System(eqs, t)

$$
\begin{align}
0 &= 1 - \left( x\left( t \right) \right)^{2} - \left( y\left( t \right) \right)^{2} \\
\frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} &= \mathtt{yˍt}\left( t \right) \\
0 &= - 2 \mathtt{xˍt}\left( t \right) x\left( t \right) - 2 y\left( t \right) \mathtt{yˍt}\left( t \right) \\
\frac{\mathrm{d} \mathtt{yˍt}\left( t \right)}{\mathrm{d}t} &= - g + y\left( t \right) \lambda\left( t \right) \\
0 &= - 2 \mathtt{xˍtt}\left( t \right) x\left( t \right) - 2 \left( \mathtt{xˍt}\left( t \right) \right)^{2} - 2 \left( \mathtt{yˍt}\left( t \right) \right)^{2} - 2 y\left( t \right) \left( - g + y\left( t \right) \lambda\left( t \right) \right)
\end{align}
$$

Note how the second and fourth equations are differential equations, and the remaining three are algebraic. This means that the system only needs two initial conditions to initialize. It will solve for the other three variables in unknowns(pend) automatically. Previously, ModelingToolkit.jl would require us to specify y and D(y) to solve the ODEProblem, since they are the chosen differential variables. However, with initialization we can provide initial conditions for (almost) any two variables. For example, we will start our pendulum perfectly horizontal:

u0 = [x => 1, D(y) => 0, g => 1]
prob = ODEProblem(pend, u0, (0.0, 10.0))
sol = solve(prob, FBDF())
plot(sol; idxs = [x, y])

Which results in:

If we try to overspecify initial conditions:

bad_u0 = [x => 1, D(y) => 0, y => 1, g => 1]
prob=ODEProblem(pend,bad_u0,(0.0,10.0))

ModelingToolkit will provide a warning when creating the ODEProblem:

┌ Warning: Initialization system is overdetermined. 3 equations for 2 unknowns

Note that not only is there one extra initial condition, but the conditions are unsatisfiable since both x and y cannot be 1 simultaneously. If we try to solve this problem, the solution fails with the InitialFailure retcode indicating an inconsistent initialization. While this model was trivial, a large hierarchical model may inadvertently introduce spurious initial conditions via defaults that initialization would detect.

Polynomialization and Homotopy Continuation

ModelingToolkit.jl is now capable of recognizing when a nonlinear system is polynomial and can target the new HomotopyContinuation.jl wrapper in NonlinearSolveHomotopyContinuation.jl. This enables finding all roots of a polynomial system.

@variables x = 0.25 y = 0.125
eqs = [0 ~ (x^2 - 5x * y + 6y^2) / (x - 0.25)
       0 ~ (x - 0.25) * (x - 0.5)]
@named sys = System(eqs, [x, y], [])
sys = complete(sys)
prob = HomotopyContinuationProblem(sys, [])
alg = HomotopyContinuationJL{true}()
sol = solve(prob, alg)::EnsembleSolution
getproperty.(sol.u, :u)

The resultant solution reports two roots:

2-element Vector{Vector{Float64}}:
 [0.5000000000000001, 0.24999999999999994]
 [0.4999999999999999, 0.1666666666666666]

Note how the input system consists of a rational function in the first equation. ModelingToolkit.jl automatically recognizes this and tracks the fact that x = 0.25 cannot be a solution. It then transforms the equation by removing the denominator to be able to target HomotopyContinuation.jl.

We can take this one step further by phrasing our system as a polynomial of non-polynomial terms. For example,

import Nemo
a = sin(x^2 - 4x + 1)
b = cos(3log(y) + 4)
eqs = [0 ~ (a^2 - 5a * b + 6b^2) / (a - 0.25)
       0 ~ (a - 0.25) * (a - 0.5)]
@named sys = System(eqs, [x, y], [])
sys = complete(sys)
prob = HomotopyContinuationProblem(sys, [])
sol = solve(prob, alg)::EnsembleSolution
getproperty.(sol.u, :u)

This is the same polynomial, except x and y are replaced by sin(x^2 - 4x + 1) and cos(3log(y) + 4) respectively. We can still use the exact same API to target the homotopy continuation solvers because ModelingToolkit recognizes that this is a polynomial on a different basis. It solves the polynomial system for a and b above, obtaining the same roots as before, and then uses the symbolic solvers in Symbolics.jl to solve for x and y from each root for a and b respectively.

4-element Vector{Vector{Float64}}:
 [3.877125135838924, 0.4208197921273323]
 [0.12287486416107604, 0.4208197921273323]
 [3.877125135838924, 0.40903223598782096]
 [0.12287486416107596, 0.40903223598782096]

We obtain four solutions here because there are two possible values of x for each value of a.

Dynamic Optimization

ModelingToolkit.jl can now perform dynamic optimization by specifying costs and constraints on an ODE.

@variables x(..) v(..)
@variables u(..) [bounds = (-1.0, 1.0), input = true]
constr = [v(1.0) ~ 0.0]
cost = [-x(1.0)]
@named block = System(
    [D(x(t)) ~ v(t), D(v(t)) ~ u(t)], t; costs = cost, constraints = constr)
block = mtkcompile(block; inputs = [u(t)])

This syntax allows targeting BoundaryValueDiffEq.jl alongside a variety of other solvers such as Pyomo, CasADi, JuMP, and InfiniteOpt.

Unified System Type

One of the largest changes in ModelingToolkit.jl version 10 was the introduction of a single, unified System type to replace the old ODESystem, SDESystem, NonlinearSystem and so on. The key insight that motivates this is that almost all systems contain the same information. Consider an ODESystem as our template. SDESystem is just an ODESystem with noise. A JumpSystem is obtained by adding jumps to an ODE, SDE or Discrete system of equations. A NonlinearSystem is structurally identical to an ODESystem with the exception of not having an independent variable. And this allows us to express more combinations of problems. For example, a dynamic optimization problem is a constrained ODE with a cost function - essentially a combination of the old ODESystem with OptimizationSystem. The new System type just validates that it has all the pieces of an ODE and additional costs/constraints. This allows growing the variety of numerical problems we can represent without a corresponding combinatorial explosion in types.

ImplicitDiscreteProblem and Improved Symbolic Callbacks

An ImplicitDiscreteProblem is the DAE analogue of DiscreteProblem where some variables can be explicitly stepped and others require solving a nonlinear system of equations. Consider the following toy example:

@variables x(t) y(t) [guess = 1.0]
k = ShiftIndex(t)
eqs = [x(k) ~ x(k-1) + x(k-2)
       exp(y) ~ x]
@mtkcompile sys = System(eqs, t)

Here x just follows the Fibonacci sequence. However, y cannot be written explicitly in terms of x by linearly rearranging the second equation. Thus, it is an implicit discrete system. ModelingToolkit.jl now supports solving such systems:

u0 = [x(k-1) => 1.0, x => 1.0]
prob = ImplicitDiscreteProblem(sys, u0, (0, 5))
sol = solve(prob, IDSolve())
sol[[x, x - exp(y)]]

Which results in:

6-element Vector{Vector{Float64}}:
 [1.0, 0.0]
 [2.0, 0.0]
 [3.0, -4.440892098500626e-16]
 [5.0, -8.881784197001252e-16]
 [8.0, -1.7763568394002505e-15]
 [13.0, 0.0]

We can see that x follows the Fibonacci sequence as expected, and the residual for the second equation is zero (within tolerance).

Why is this functionality useful? This now allows us to represent symbolic callbacks in a declarative manner. Going back to our pendulum:

@parameters g
@variables x(t)
@variables y(t) [state_priority = 10, guess = 1.0]
@variables λ(t) [guess = 1.0]
eqs = [D(D(x)) ~ λ * x
       D(D(y)) ~ λ * y - g
       x^2 + y^2 ~ 1]
devts = [2 => [x ~ Pre(x) - 0.1 * sign(Pre(x))]]
@mtkcompile pend = System(eqs, t;
    discrete_events = devts)

Here we’ve added an event that triggers every two seconds, and decreases the value of x by 10%. With ModelingToolkit.jl version 10, Pre is used to refer to the value of a variable before the callback. Note that in the equations of the pendulum, x is an algebraic variable. If we were to implement this without ModelingToolkit.jl and simply change the value of x in a callback, the DAE solver would detect that the algebraic equations are not satisfied. Since the solver is unaware that we intentionally changed x, it will simply fix the values of the differential variables and solve for the algebraic ones. Thus, the value of x will be reset to the manifold. This is a common issue that makes it seem as if an event did not trigger when in fact the DAE solver simply tried its best to solve the problem. ModelingToolkit.jl phrases the update as an implicit discrete system by including the given update equation alongside the algebraic equations of the system. Since Pre(x) is fixed, the value of x is fixed and the callback solves for the rest of the variables to re-satisfy the DAE constraints. If we create and solve the ODEProblem as below:

u0 = [x => -1/2, D(x) => 1/2, g => 1]
prob = ODEProblem(pend, u0, (0.0, 10.0))
sol = solve(prob, FBDF())

We can see that the values of x and y change suddenly every 2s.

And more

There are many features that are not covered here. Check out the talk or ModelingToolkit.jl documentation for more information.

Future Goals

• Analytical integration

  We aim to enable symbolically integrating simple differential equations. This will both improve solve times as well as provide more accurate results. For example, the equation D(x) ~ a * x can be automatically integrated to x ~ x0 * exp(a * t).

• Symbolic structs

  Support for proper symbolic structs would allow more concise representation of information and the use of getproperty to access fields, similar to normal Julia structs.

• Improved support for Array equations

  While ModelingToolkit.jl is capable of handling arrays in some forms, this functionality can be extended. Notably, we aim to allow simplification of systems with array equations without requiring scalarization of the equations. This can significantly improve simplification times for large systems.

• Integration with Reactant.jl

  With support for array equations, we aim to target optimizing compilers such as Reactant.jl to be able to generate highly efficient code for array operations.

• Rewrite of SymbolicUtils.jl

  A comprehensive overhaul of the underlying symbolic infrastructure which will significantly improve performance of symbolic operations is almost complete.

• Specialized code generation to exploit common system structures

  We aim to automatically detect and generate code for common structures in equations. For example, many equations can be phrased as a linear transformation of the state with some additional nonlinearities. Typically, the linear part of the system is stiff and the nonlinearities are non-stiff. The work-in-progress SemilinearODEProblem aims to recognize such problems of the form D(x)[1:n] ~ A[1:n, 1:n] * x[1:n] + C(x)[1:n] and target IMEX solvers in the OrdinaryDiffEq ecosystem, enabling faster simulation.