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:

using ModelingToolkit, NonlinearSolve

@variables x y z
@parameters σ ρ β

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

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

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

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

retcode: Success
u: 3-element Vector{Float64}:
  0.0
 -5.861075866041436e-17
 -2.1979034497655384e-17

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:

equations(ns)

\[ \begin{align} 0 =& \left( - x + y \right) \sigma \\ 0 =& - y + x \left( - z + \rho \right) \\ 0 =& x y - z \beta \end{align} \]

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

calculate_jacobian(ns)

\[ \begin{equation} \left[ \begin{array}{ccc} - \sigma & \sigma & 0 \\ - z + \rho & -1 & - x \\ y & x & - \beta \\ \end{array} \right] \end{equation} \]

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

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

retcode: Success
u: 3-element Vector{Float64}:
  0.0
 -5.861075866041436e-17
 -2.1979034497655384e-17

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:

@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),
]
@named sys = NonlinearSystem(eqs, [u1, u2, u3, u4, u5], [])

\[ \begin{align} 0 =& u1 - \sin\left( u5 \right) \\ 0 =& u2 - \cos\left( u1 \right) \\ 0 =& u3 - \mathrm{hypot}\left( u1, u2 \right) \\ 0 =& u4 - \mathrm{hypot}\left( u2, u3 \right) \\ 0 =& u5 - \mathrm{hypot}\left( u4, u1 \right) \end{align} \]

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

sys = structural_simplify(sys)

\[ \begin{align} 0 =& u5 - \mathrm{hypot}\left( u4, u1 \right) \end{align} \]

equations(sys)

\[ \begin{align} 0 =& u5 - \mathrm{hypot}\left( u4, u1 \right) \end{align} \]

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

observed(sys)

\[ \begin{align} u1 =& \sin\left( u5 \right) \\ u2 =& \cos\left( u1 \right) \\ u3 =& \mathrm{hypot}\left( u1, u2 \right) \\ u4 =& \mathrm{hypot}\left( u2, u3 \right) \end{align} \]

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

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

retcode: Success
u: 1-element Vector{Float64}:
 1.6069926947050053

We can then use symbolic indexing to retrieve any variable:

sol[u1]

0.9993449829954304

sol[u2]

0.540853367725015

sol[u3]

1.1363154317431527

sol[u4]

1.2584653852200776

sol[u5]

1.6069926947050053

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.