# ModelingToolkit vs plain DifferentialEquations

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322>\
**Category:** Modelling & Simulations\
**Tags:** differentialequation\
**Created:** [May 19, 2022, 6:26pm UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322 "2022-05-19T18:26:16Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![mzhenirovskyy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzhenirovskyy/32/7549_2.png) [@mzhenirovskyy](https://discourse.julialang.org/u/mzhenirovskyy)\
**Post date:** [May 19, 2022, 6:26pm UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/1 "2022-05-19T18:26:16Z")

</div>

Hi All,

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

When I have the DAE system in the next form:

```julia
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:

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

retcode: InitialFailure

```

But when I use the ModelingToolkit:

```julia
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.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [May 19, 2022, 11:09pm UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/2 "2022-05-19T23:09:02Z")

</div>

Aren’t you missing the derivative terms in the DifferentialEquations.jl version? See here for an example: [Differential Algebraic Equations · DifferentialEquations.jl](https://diffeq.sciml.ai/latest/tutorials/dae_example/#Implicitly-Defined-Differential-Algebraic-Equations-(DAEs))

---

<div class="post-metadata">

**Author:** ![mzhenirovskyy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzhenirovskyy/32/7549_2.png) [@mzhenirovskyy](https://discourse.julialang.org/u/mzhenirovskyy)\
**Post date:** [May 19, 2022, 11:58pm UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/3 "2022-05-19T23:58:01Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [May 20, 2022, 12:10am UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/4 "2022-05-20T00:10:18Z")

</div>

> [@mzhenirovskyy](#):
>
> `out[1] = u[8] / 1.072`

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

```julia
 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

```

---

<div class="post-metadata">

**Author:** ![mzhenirovskyy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzhenirovskyy/32/7549_2.png) [@mzhenirovskyy](https://discourse.julialang.org/u/mzhenirovskyy)\
**Post date:** [May 20, 2022, 12:39am UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/5 "2022-05-20T00:39:54Z")

</div>

Thank you @isaacsas!  
I was stupid…

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [May 20, 2022, 12:43am UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/6 "2022-05-20T00:43:09Z")

</div>

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](https://github.com/SciML/DiffEqDocs.jl/issues)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 20, 2022, 2:00am UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/7 "2022-05-20T02:00:45Z")

</div>

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](https://mtk.sciml.ai/v8.11/mtkitize_tutorials/modelingtoolkitize_index_reduction/), 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.

---

<div class="post-metadata">

**Author:** ![mzhenirovskyy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzhenirovskyy/32/7549_2.png) [@mzhenirovskyy](https://discourse.julialang.org/u/mzhenirovskyy)\
**Post date:** [May 20, 2022, 2:15am UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/8 "2022-05-20T02:15:05Z")

</div>

@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.

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 20, 2022, 2:19am UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/9 "2022-05-20T02:19:28Z")

</div>

> [@mzhenirovskyy](#):
>
> Can I use ModelingToolkit together with DiffEqFlux? - my final goal is to run optimization based on DAE. Do you have an example?

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).

> [@mzhenirovskyy](#):
>
> 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.

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?

---

<div class="post-metadata">

**Author:** ![mzhenirovskyy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzhenirovskyy/32/7549_2.png) [@mzhenirovskyy](https://discourse.julialang.org/u/mzhenirovskyy)\
**Post date:** [May 20, 2022, 2:33am UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/10 "2022-05-20T02:33:04Z")

</div>

@ChrisRackauckas IDA says

```julia
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})

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 20, 2022, 2:58am UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/11 "2022-05-20T02:58:51Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 20, 2022, 3:33am UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/12 "2022-05-20T03:33:31Z")

</div>

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

```julia
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](https://github.com/SciML/DiffEqBase.jl/pull/752)

---

<div class="post-metadata">

**Author:** ![mzhenirovskyy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzhenirovskyy/32/7549_2.png) [@mzhenirovskyy](https://discourse.julialang.org/u/mzhenirovskyy)\
**Post date:** [May 20, 2022, 3:42am UTC](https://discourse.julialang.org/t/modelingtoolkit-vs-plain-differentialequations/81322/13 "2022-05-20T03:42:45Z")

</div>

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