Symbolic Nonlinear System Definition and Acceleration via ModelingToolkit

ModelingToolkit.jl is a symbolic-numeric modeling system for the Julia SciML ecosystem. It adds a high-level interactive interface for the numerical solvers which can make it easy to symbolically modify and generate equations to be solved. The basic form of using ModelingToolkit looks as follows:

import ModelingToolkit as MTK
import NonlinearSolve as NLS

MTK.@variables x y z
MTK.@parameters σ ρ β

# Define a nonlinear system
eqs = [0 ~ σ * (y - x), 0 ~ x * (ρ - z) - y, 0 ~ x * y - β * z]
MTK.@mtkbuild ns = MTK.NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])

u0 = [x => 1.0, y => 0.0, z => 0.0]

ps = [σ => 10.0
      ρ => 26.0
      β => 8 / 3]

prob = NLS.NonlinearProblem(ns, u0, ps)
sol = NLS.solve(prob, NLS.NewtonRaphson())

Symbolic Derivations of Extra Functions

As a symbolic system, ModelingToolkit can be used to represent the equations and derive new forms. For example, let's look at the equations:

MTK.equations(ns)

We can ask it what the Jacobian of our system is via calculate_jacobian:

MTK.calculate_jacobian(ns)

We can tell MTK to generate a computable form of this analytical Jacobian via jac = true to help the solver use efficient forms:

prob = NLS.NonlinearProblem(ns, u0, ps, jac = true)
sol = NLS.solve(prob, NLS.NewtonRaphson())

Symbolic Simplification of Nonlinear Systems via Tearing

One of the major reasons for using ModelingToolkit is to allow structural simplification of the systems. It's very easy to write down a mathematical model that, in theory, could be solved more simply. Let's take a look at a quick system:

MTK.@variables u1 u2 u3 u4 u5
eqs = [0 ~ u1 - sin(u5), 0 ~ u2 - cos(u1), 0 ~ u3 - hypot(u1, u2),
    0 ~ u4 - hypot(u2, u3), 0 ~ u5 - hypot(u4, u1)]
MTK.@named sys = MTK.NonlinearSystem(eqs, [u1, u2, u3, u4, u5], [])

If we run structural simplification, we receive the following form:

sys = MTK.structural_simplify(sys)

MTK.equations(sys)

How did it do this? Let's look at the observed to see the relationships that it found:

MTK.observed(sys)

Using ModelingToolkit, we can build and solve the simplified system:

u0 = [u5 .=> 1.0]
prob = NLS.NonlinearProblem(sys, u0)
sol = NLS.solve(prob, NLS.NewtonRaphson())

We can then use symbolic indexing to retrieve any variable:

sol[u1]

sol[u2]

sol[u3]

sol[u4]

sol[u5]

Component-Based and Acausal Modeling

If you're interested in building models in a component or block based form, such as seen in systems like Simulink or Modelica, take a deeper look at ModelingToolkit.jl's documentation which goes into detail on such features.