ModelingToolkit IR

ModelingToolkit IR, which falls under the Expression abstract type, mirrors the Julia AST but allows for easy mathematical manipulation by itself following mathematical semantics. The base of the IR is the Variable type, which defines a symbolic variable. These variables are combined using Operations, which are registered functions applied to the various variables. These Operations then perform automatic tracing, so normal mathematical functions applied to an Operation generate a new Operation. For example, op1 = x+y is one Operation and op2 = 2z is another, and so op1*op2 is another Operation. Then, at the top, an Equation, normally written as op1 ~ op2, defines the symbolic equality between two operations.


struct Variable{T} <: Function

A named variable which represents a numerical value. The variable is uniquely identified by its name, and all variables with the same name are treated as equal.


  • name

    The variable's unique name.

For example, the following code defines an independent variable t, a parameter α, a function parameter σ, a variable x, which depends on t, a variable y with no dependents, a variable z, which depends on t, α, and x(t) and parameters β₁ and β₂.

t = Variable(:t)()  # independent variables are treated as known
α = Variable(:α)()  # parameters are known
σ = Variable(:σ)    # left uncalled, since it is used as a function
w = Variable(:w)   # unknown, left uncalled
x = Variable(:x)(t)  # unknown, depends on `t`
y = Variable(:y)()   # unknown, no dependents
z = Variable(:z)(t, α, x)  # unknown, multiple arguments
β₁ = Variable(:β, 1)() # with index 1
β₂ = Variable(:β, 2)() # with index 2

expr = β₁ * x + y^α + σ(3) * (z - t) - β₂ * w(t - 1)
struct Operation <: Expression

An expression representing the application of a function to symbolic arguments.


  • op

    The function to be applied.

  • args

    The arguments the function is applied to.


Operations can be built by application of most built-in mathematical functions to other Expression instances:

julia> using ModelingToolkit

julia> @variables x y;

julia> op1 = sin(x)

julia> typeof(op1.op)

julia> op1.args
1-element Array{Expression,1}:

julia> op2 = x + y
x() + y()

julia> typeof(op2.op)

julia> op2.args
2-element Array{Expression,1}:
struct Equation

An equality relationship between two expressions.


  • lhs

    The expression on the left-hand side of the equation.

  • rhs

    The expression on the right-hand side of the equation.


Function Registration

The ModelingToolkit graph only allowed for registered Julia functions for the operations. All other functions are automatically traced down to registered functions. By default, ModelingToolkit.jl pre-registers the common functions utilized in the AD package ruleset DiffRules.jl and pre-defines their derivatives. However, the user can utilize the @register macro to add their function to allowed functions of the computation graph.


Registers a function call as a primitive for the Operation graph of the ModelingToolkit IR. Example:

@register f(x,y)

registers f as a possible two-argument function.

You may also want to tell ModelingToolkit the derivative of the registered function. Here is an example to do it

julia> using ModelingToolkit

julia> foo(x, y) = sin(x) * cos(y)
foo (generic function with 1 method)

julia> @parameters t; @variables x(t) y(t) z(t); @derivatives D'~t;

julia> @register foo(x, y)
foo (generic function with 4 methods)

julia> foo(x, y)
foo(x(t), y(t))

julia> ModelingToolkit.derivative(::typeof(foo), (x, y), ::Val{1}) = cos(x) * cos(y) # derivative w.r.t. the first argument

julia> ModelingToolkit.derivative(::typeof(foo), (x, y), ::Val{2}) = -sin(x) * sin(y) # derivative w.r.t. the second argument

julia> isequal(expand_derivatives(D(foo(x, y))), expand_derivatives(D(sin(x) * cos(y))))

Derivatives and Differentials

A Differential(op) is a partial derivative with respect to the operation op, which can then be applied to some other operations. For example, D=Differential(t) is what would commonly be referred to as d/dt, which can then be applied to other operations using its function call, so D(x+y) is d(x+y)/dt.

By default, the derivatives are left unexpanded to capture the symbolic representation of the differential equation. If the user would like to expand out all of the differentials, the expand_derivatives function eliminates all of the differentials down to basic one-variable expressions.

derivative(O, idx)

Calculate the derivative of the op O with respect to its argument with index idx.


julia> using ModelingToolkit

julia> @variables x y;

julia> ModelingToolkit.derivative(sin(x), 1)

Note that the function does not recurse into the operation's arguments, i.e., the chain rule is not applied:

julia> myop = sin(x) * y^2
sin(x()) * y() ^ 2

julia> typeof(myop.op)  # Op is multiplication function

julia> ModelingToolkit.derivative(myop, 1)  # wrt. sin(x)
y() ^ 2

julia> ModelingToolkit.derivative(myop, 2)  # wrt. y^2
struct Differential <: Function

Represents a differential operator.


  • x

    The variable or expression to differentiate with respect to.


julia> using ModelingToolkit

julia> @variables x y;

julia> D = Differential(x)

julia> D(y)  # Differentiate y wrt. x

Note that the generation of sparse matrices simply follows from the Julia semantics imbued on the IR, so sparse(jac) changes a dense Jacobian to a sparse Jacobian matrix.

Adding Derivatives

There is a large amount of derivatives pre-defined by DiffRules.jl. Note that Expression types are defined as <:Real, and thus any functions which allow the use of real numbers can automatically be traced by the derivative mechanism. Thus, for example:

f(x,y,z) = x^2 + sin(x+y) - z

automatically has the derivatives defined via the tracing mechanism. It will do this by directly building the operation the internals of your function and differentiating that.

However, in many cases you may want to define your own derivatives so that way automatic Jacobian etc. calculations can utilize this information. This can allow for more succinct versions of the derivatives to be calculated in order to better scale to larger systems. You can define derivatives for your own function via the dispatch:

# `N` arguments are accepted by the relevant method of `my_function`
ModelingToolkit.derivative(::typeof(my_function), args::NTuple{N,Any}, ::Val{i})

where i means that it's the derivative of the ith argument. args is the array of arguments, so, for example, if your function is f(x,t), then args = [x,t]. You should return an Operation for the derivative of your function.

For example, sin(t)'s derivative (by t) is given by the following:

ModelingToolkit.derivative(::typeof(sin), args::NTuple{1,Any}, ::Val{1}) = cos(args[1])

IR Manipulation

ModelingToolkit.jl provides functionality for easily manipulating Expression types. Most of the functionality comes by the Expression type obeying the standard mathematical semantics. For example, if one has A a matrix of Expression, then A^2 calculates the Expressions for the squared matrix. In that sense, it is encouraged that one uses standard Julia for performing a lot of the manipulation on the IR, as, for example, calculating the sparse form of the matrix via sparse(A) is valid, legible, and easily understandable to all Julia programmers.

Other additional manipulation functions are given below.

Missing docstring.

Missing docstring for simplify_constants. Check Documenter's build log for details.

Missing docstring.

Missing docstring for get_variables. Check Documenter's build log for details.

Missing docstring.

Missing docstring for substitute_expr!. Check Documenter's build log for details.

Expression Generation and build_function

At any time, Julia expressions can be generated from ModelingToolkit IR by using convert(Expr,x). This performs some cleaning to return an expression without extraneous pieces that commonly matches expressions one would write in functions like those for differential equation solvers and optimization libraries.

Additionally, the core compilation process of ModelingToolkit IR is build_function. build_function takes an operation or an AbstractArray of operations and generates a compilable version of the model for numerical solvers.



Generates a numerically-usable function from a ModelingToolkit Expression. If the Expression is an Operation, the generated function is a function with a scalar output, otherwise if it's an AbstractArray{Operation}, the output is two functions, one for out-of-place AbstractArray output and a second which is a mutating function. The outputted functions match the given argument order, i.e., f(u,p,args...) for the out-of-place and scalar functions and f!(du,u,p,args..) for the in-place version.

build_function(ex, args...;
               conv = simplified_expr, expression = Val{true},
               checkbounds = false, convert_oop = true,
			   force_SA = false,
               linenumbers = false, target = JuliaTarget())


  • ex: The Expression to compile
  • vs: The variables of the expression
  • ps: The parameters of the expression
  • args: Extra arguments to the function
  • conv: The conversion function of the Operation to Expr. By default this uses the simplified_expr function utilized in convert(Expr,x).
  • expression: Whether to generate code or whether to generate the compiled form. By default, expression = Val{true}, which means that the code for the function is returned. If Val{false}, then the returned value is a compiled Julia function, which utilizes GeneralizedGenerated.jl in order to world-age free.

Keyword Arguments:

  • checkbounds: For whether to enable bounds checking inside of the generated function. Defaults to false, meaning that @inbounds is applied.
  • linenumbers: Determines whether the generated function expression retains the line numbers. Defaults to true.
  • convert_oop: Determines whether the OOP version should try to convert the output to match the type of the first input. This is useful for cases like LabelledArrays or other array types that carry extra information. Defaults to true.
  • force_SA: Forces the output of the OOP version to be a StaticArray. Defaults to false, and outputs a static array when the first argument is a static array.
  • target: The output target of the compilation process. Possible options are:
    • JuliaTarget: Generates a Julia function
    • CTarget: Generates a C function
    • StanTarget: Generates a function for compiling with the Stan probabilistic programming language
    • MATLABTarget: Generates an anonymous function for use in MATLAB and Octave environments