ModelingToolkit vs plain DifferentialEquations

Specific Domains Modelling & Simulations
1

Hi All,

What is the magic that happens inside of the “ModelingToolkit.ODESystem”?

When I have the DAE system in the next form:

using DifferentialEquations, OrdinaryDiffEq, Sundials


u0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
du0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

t1 = [0.0:0.1:10.0;]
tspan = (t1[1], last(t1))

function dfdt(out, du, u, p, t)
    out[1] = u[8] / 1.072
    out[2] = u[9] / 0.2638
    out[3] = u[10] / 1.682
    out[4] = u[11] / 0.8365
    out[5] = u[12] / 0.7265
    out[6] = u[2] / 1.304
    out[7] = u[4] / 0.8586
    out[8] = u[10] + u[11] + u[7] - u[9] - u[6]
    out[9] = u[1] - u[3] - u[2]
    out[10] = u[3] - u[5] - u[4]
    out[11] = u[8] + u[9] + u[6] - 1.0 + u[1]
    out[12] = u[11] + u[7] - u[12] - u[5]
end

differential_vars = [true, true, true, true, true, true, true,
    false, false, false, false, false]

prob = DAEProblem(dfdt, du0, u0, tspan, differential_vars=differential_vars)
sol = solve(prob, IDA())

I can’t get a reasonable solution.
The result:

[IDAS ERROR]  IDACalcIC
  Newton/Linesearch algorithm failed to converge.

retcode: InitialFailure

But when I use the ModelingToolkit:

using ModelingToolkit, DifferentialEquations,Sundials


@variables t x[1:12](t)
D = Differential(t)

eqs = [
      D(x[1]) ~ x[8] / 1.072,
      D(x[2]) ~ x[9] / 0.2638,
      D(x[3]) ~ x[10] / 1.682,
      D(x[4]) ~ x[11] / 0.8365,
      D(x[5]) ~ x[12] / 0.7265,
      D(x[6]) ~ x[2] / 1.304,
      D(x[7]) ~ x[4] / 0.8586,
      0 ~ x[10] + x[11] + x[7] - x[9] - x[6],
      0 ~ x[1] - x[3] - x[2],
      0 ~ x[3] - x[5] - x[4],
      0 ~ x[8] + x[9] + x[6] - 1.0 + x[1],
      0 ~ x[11] + x[7] - x[12] - x[5]
]

@named sys = ODESystem(eqs)

u0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
du0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

t1 = [0.0:0.1:10.0;]
tspan = (t1[1], last(t1))

prob = DAEProblem(sys, du0, u0, tspan)
sol = solve(prob, IDA())

I can get a reasonable result.

The solver (IDA()) and DAEs are absolutely the same.

2

Aren’t you missing the derivative terms in the DifferentialEquations.jl version? See here for an example: Differential Algebraic Equations · DifferentialEquations.jl

1 Like
3

@isaacsas
I didn’t get you. Could you point me to the derivative terms you mean?

4

Take a look at that example I linked, it shows how to write DAEs in DifferentialEquations.jl. Your dfdt function is just algebraic equations currently, it isn’t using the derivatives, du. This is what you should be using

 function dfdt(out, du, u, p, t)
           out[1] = u[8] / 1.072 - du[1]
           out[2] = u[9] / 0.2638 - du[2]
           out[3] = u[10] / 1.682 - du[3]
           out[4] = u[11] / 0.8365 - du[4]
           out[5] = u[12] / 0.7265 - du[5]
           out[6] = u[2] / 1.304 - du[6]
           out[7] = u[4] / 0.8586 - du[7]
           out[8] = u[10] + u[11] + u[7] - u[9] - u[6]
           out[9] = u[1] - u[3] - u[2]
           out[10] = u[3] - u[5] - u[4]
           out[11] = u[8] + u[9] + u[6] - 1.0 + u[1]
           out[12] = u[11] + u[7] - u[12] - u[5]
       end
5

Thank you @isaacsas!
I was stupid…

6

No worries! Glad we got it fixed!

With that said, as someone new to using the DAE stuff if anything in the docs isn’t clear please let us know / make suggestions! You can open suggestions on ways to improve the docs at https://github.com/SciML/DiffEqDocs.jl/issues

7

To answer the other bits of the question though, DAEs are hard and can be very sensitive to the exact formulation of the DAE that you give. For example, out[9] = u[1] - u[3] - u[2] defines an algebraic relationship (0 = u[1] - u[3] - u[2]), but one could have also done u[2] = u[1] - u[3], substituted out u[2], and received a smaller set of equations. In many cases, this can improve the stability of the solving by enforcing the equality exactly and structurally instead of approximately. Removing variables can also of course make the computation faster by making the Jacobians smaller.

However, there are many details in just that: why remove u[2] instead of u[3]? This is the question of state selection. Some choices are stable, while other choices might not be. But there are other effects as well. DAEs have something known as a differential index, and DAE solvers require this differential index to be 1 (well, it can be 2 sometimes, maybe 3, details details). You can symbolically transform equations via the Pantelides algorithm, choosing some equations to differentiate and then substitute in values to receive a mathematically-equivalent system that is numerically more robust. This is seen for example on a Cartesian pendulum example in the documentation, where the original formulation is unstable but the transformed equations are stable. Then there’s the question of tearing nonlinear equations for stability and performance.

When you dig deep down into it, you’ll see that what emerges is that structural_simplify puts all of these concepts together to generate more efficient and stable numerical formulations of DAEs, effectively making the job of the DAE solver a lot easier. For this reason, it’s highly recommended for “most” DAE solves. That said, at the end of the day it’s just a symbolic front-end for DifferentialEquations.jl, so you can always write the code manually, and right now it has some limitations, notably the code it generates always scalarizes loops and there’s limited support for event handling (both are being actively worked on, but are fairly difficult issues). Thus for the most part, we suggest that people doing DAE modeling should be using ModelingToolkit, especially because it opens up new problem formulations such as the ODAEProblem which requires nonlinear tearing to even describe the code it generates (i.e. it’s pretty much a symbolic-only variation of the problem), and there’s the soon to come SemilinearODEProblem for accelerating many situations as well using generate SplitODEProblem code. But if you want the bare metal, use DifferentialEquations.jl since that’s what it is always targeting anyways.

2 Likes
8

@isaacsas, You are right the most significant disadvantage of Julia is their lack of good documentation. Direct competitors like Python/Matlab have pretty good docs.

@ChrisRackauckas, Yes, I’m going to use ModelingToolkit. Can I use ModelingToolkit together with DiffEqFlux? - my final goal is to run optimization based on DAE. Do you have an example?

Can you guys point me in the right direction on fixing the solver’s “Instability detected”?

I continue to play with the same system, but now in the form of ODE with the Mass Matrix.

using DifferentialEquations, OrdinaryDiffEq
using Sundials, DASSL, DASKR

function dfdt!(du, u, p, t)
      du[1] = u[8] / 1.072
      du[2] = u[9] / 0.2638
      du[3] = u[10] / 1.682
      du[4] = u[11] / 0.8365
      du[5] = u[12] / 0.7265
      du[6] = u[2] / 1.304
      du[7] = u[4] / 0.8586
      du[8] = u[10] + u[11] + u[7] - u[9] - u[6]
      du[9] = u[1] - u[3] - u[2]
      du[10] = u[3] - u[5] - u[4]
      du[11] = u[8] + u[9] + u[6] - 1.0 + u[1]
      du[12] = u[11] + u[7] - u[12] - u[5]
      return nothing
end

u0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
M = Diagonal([1.,1.,1.,1.,1.,1.,1.,0.,0.,0.,0.,0.])

t1 = [0.0:0.1:10.0;]
tspan = (t1[1], last(t1))

f = ODEFunction(dfdt!, mass_matrix=M)
prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob, Rodas5())

The result is “Instability detected. Aborting”
I tried the different solvers: Rodas5, Rosenbrock23, RadauIIA5, ROS3P, Rodas3, RosShamp4, Veldd4, Velds4, QBDF, dassl, daskr - the same result.

9

Generally the easiest way to do this is to just not have MTK in the optimization loop. Generate the equations to solve once, and then use remake to keep updating the parameters in the optimization loop. Since it’s not even in the optimization part, it’s the same as any other DifferentialEquations.jl code (and avoids any recompilation).

This is likely one of the stability issues I mentioned. By eye, the DAE does not look like it’s index-2, but just having a bunch of linear equality relationships will decrease the stability and make it much easier to have a singular Jacobian pop up, which then would give instability and abort. I’d have to dig in more to see, but starting at all zeros on that equation does not give me confidence of a nonsingular Jacobian. What did IDA throw?

10

@ChrisRackauckas IDA says

ERROR: MethodError: no method matching __init(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(dfdt!), Diagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::IDA{:Dense, Nothing, Nothing}, ::Vector{Any}, ::Vector{Any}, ::Vector{Any})
11

IDA needs a DAEProblem. Oh I see what you mean, this mass matrix formulation is doing something a little different. I can’t check it right now, but see if the simplified MTK generated ODEProblem solves it. If so, it may be due to the initialization that it’s finding. Also try initializationalg = ShampineCollocation() as a keyword argument.

1 Like
12

Note I went ahead and cleaned up that error message. It should read now like:

julia> sol = solve(prob, IDA())
ERROR: Incompatible problem+solver pairing.
For example, this can occur if an ODE solver is passed with an SDEProblem.
Solvers are only capable of handling specific problem types. Please double
check that the chosen pairing is capable for handling the given problems.

Problem type: ODEProblem
Solver type: IDA
Problem types compatible with the chosen solver: DAEProblem

Stacktrace:
 [1] check_prob_alg_pairing(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(dfdt!), Diagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::IDA{:Dense, Nothing, Nothing})
   @ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:505
 [2] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(dfdt!), Diagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::SciMLBase.NullParameters, args::IDA{:Dense, Nothing, Nothing}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:270
 [3] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(dfdt!), Diagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::SciMLBase.NullParameters, args::IDA{:Dense, Nothing, Nothing})
   @ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:260
 [4] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(dfdt!), Diagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::IDA{:Dense, Nothing, Nothing}; sensealg::Nothing, u0::Nothing, p::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:255
 [5] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(dfdt!), Diagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::IDA{:Dense, Nothing, Nothing})
   @ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:250
 [6] top-level scope
   @ REPL[1]:1

after https://github.com/SciML/DiffEqBase.jl/pull/752

1 Like
13

@ChrisRackauckas initializationalg = ShampineCollocationInit() doesn’t help, the same output.
I will try to play with MTK and the initial conditions.