ModelingToolkit.jl and Dyad

All code in the Dyad language lowers to Julia code that uses ModelingToolkit.jl. If you inspect the generated folder in your Dyad project, for example, you will see all your components and analyses lowered to ModelingToolkit.jl code.

Because Dyad is transparently lowered to ModelingToolkit, you can use all the power of ModelingToolkit, the SciML suite of tools, and the Julia ecosystem to programmatically interact with your models.

Let's start with an introduction to ModelingToolkit and how it works. Then, we'll show how use ModelingToolkit to programmatically interact with your Dyad models.

An introduction to ModelingToolkit

The ModelingToolkit ecosystem is built on top of the DifferentialEquations.jl suite of solvers, and the Symbolics.jl suite of symbolic tools.

In order to run this tutorial, you'll need the following packages:

using ModelingToolkit,
      Symbolics,              # Symbolic expressions and manipulation
      ForwardDiff,            # Compute the derivative of a function
      OrdinaryDiffEqDefault   # Solve ODEs

Symbolic expressions

ModelingToolkit is based on Symbolics.jl. The main unit of Symbolics.jl is the Variable, which you can define by the @variables macro:

@variables t

$$\left[\begin{array}{c}t\end{array}\right]$$

A Symbolics variable can also depend on an independent variable, so you can say that x is a function of t:

@variables x(t)

$$\left[\begin{array}{c}x\left(t\right)\end{array}\right]$$

This can also go multidimensional, so you can say that z is a function of x, y, and t:

@variables x1 y1 t1
@variables z(x1, y1, t1)

$$\left[\begin{array}{c}z(\mathtt{x}\mathtt{1},\mathtt{y}\mathtt{1},\mathtt{t}\mathtt{1})\end{array}\right]$$

In Symbolics, you define equations by the ~ operator (not =, since that becomes a Julia variable assignment). For example, let's define an equation for sin(t):

eq = x ~ sin(t)

$$x\left(t\right)=\sin \left(t\right)$$

These equations can be as complex as you want, and can even include function calls provided that Symbolics.jl can trace through them. We will explain more about tracing two sections from now.

Taking derivatives

Symbolics also allows you to define operators. For example, you can define an operator that is the differential of a variable:

D = Differential(t)

Differential(t)

If we want to take the derivative of sin(t) with respect to t, we can call the differential operator on sin(t).

eq_prime = D(sin(t))

$$\frac{\mathrm{d}\sin \left(t\right)}{\mathrm{d}t}$$

This is itself a symbolic expression - you do not get the result immediately. But you can expand it and get the result:

expand_derivatives(eq_prime)

$$\cos \left(t\right)$$

This is what we expect - cos(t).

Equations to numerical functions

Now, let's go into the interesting bit - generating fast numerical functions from our symbolic expressions.

Let's create some functions that can be solved numerically. We'll start with a simple one:

@variables y(t)

eqs = [
    sin(2x)
    x + 2y
]

$$\left[\begin{array}{c}\sin \left(2x\left(t\right)\right)\\ x\left(t\right)+2y\left(t\right)\end{array}\right]$$

We pass the array of expressions/equations (they do not have to be equations) to Symbolics.build_function, which will give us an out-of-place (op, returns an array) function and an in-place (ip, mutates the input) function.

f_op, f_ip = build_function(eqs, [x, y]; expression = false) # expression=false returns a RuntimeGeneratedFunction
f_op

RuntimeGeneratedFunction(#=in Symbolics=#, #=using Symbolics=#, :((ˍ₋arg1,)->#= /home/actions-runner-10/.julia/packages/Symbolics/JCopU/src/build_function.jl:366 =# @inbounds(begin
              #= /home/actions-runner-10/.julia/packages/Symbolics/JCopU/src/build_function.jl:366 =#
              begin
                  #= /home/actions-runner-10/.julia/packages/SymbolicUtils/0GKgW/src/code.jl:409 =#
                  #= /home/actions-runner-10/.julia/packages/SymbolicUtils/0GKgW/src/code.jl:410 =#
                  #= /home/actions-runner-10/.julia/packages/SymbolicUtils/0GKgW/src/code.jl:411 =#
                  begin
                      begin
                          #= /home/actions-runner-10/.julia/packages/SymbolicUtils/0GKgW/src/code.jl:510 =#
                          (SymbolicUtils.Code.create_array)(typeof(ˍ₋arg1), nothing, Val{1}(), Val{(2,)}(), NaNMath.sin((*)(2, ˍ₋arg1[1])), (+)(ˍ₋arg1[1], (*)(2, ˍ₋arg1[2])))
                      end
                  end
              end
          end)))

Let's invoke this function:

f_op([1, 2])

2-element Vector{Float64}:
 0.9092974268256817
 5.0

We can also invoke the in-place function on an array, whose values it will mutate:

du = [0.0, 0.0]
f_ip(du)
du

[0.0, 0.0]

Other cool things you can do with Symbolics.jl:

simplify(eqs)
substitute(eqs, [x => t])
symbolic_solve

Tracing and registered functions

In general, Symbolics can trace through most simple mathematical Julia expressions, like the function below.

function g(x, y)
    return 2x + 3y
end

g (generic function with 1 method)

You can call this with numbers, of course,

g(1, 2)

8

and if you call it with Symbolics variables, it will return a Symbolics expression.

g(x, y)

$$2x\left(t\right)+3y\left(t\right)$$

If your function uses your input in some non-traceable or non-mathematical way, you will have to "register" it, which tells Symbolics to treat it as a black box and not try to trace through it.

# models N*x^2
function h(x, N)
    z = 0
    for i=1:N
        z += x
    end
    return z*x
end

h (generic function with 1 method)

h(2, 2)

8

This will error, because you are trying to use a symbolic variable in a range.

h(x, y)

TypeError: non-boolean (Symbolics.Num) used in boolean context
A symbolic expression appeared in a Boolean context. This error arises in situations where Julia expects a Bool, like
if boolean_condition		 use ifelse(boolean_condition, then branch, else branch)
x && y				 use x & y
boolean_condition ? a : b	 use ifelse(boolean_condition, a, b)
but a symbolic expression appeared instead of a Bool. For help regarding control flow with symbolic variables, see https://docs.sciml.ai/ModelingToolkit/dev/basics/FAQ/#How-do-I-handle-if-statements-in-my-symbolic-forms?

Instead, you can inspect this:

@register_symbolic h(x, y)
h(x, y)

$$h(x\left(t\right),y\left(t\right))$$