I’m trying out ModelingToolkit.jl and want to investigate what is currently possible and what not and compare it to the OpenModelica Compiler Backend.
Let’s say we have the following DAE:
der(x) = y + a y = 3*x der(y) = a*x b = a/5
which should be simple but interesting example because it has:
- Differential index 2 (Maximum number of differentiations of all equations needed to get a solvable ODE).
- A linear loop.
- Algebraic variables.
To get an executable simulation from this OpenModelica will do:
- Matching and sorting
- Index reduction with dummy-states
See paper The OpenModelica Integrated Environment for Modeling, Simulation, and Model-Based Development chapter 4 on how the OpenModelica Backend is working.
I’m trying to solve the example from above with ModelingToolkit.jl (v6.5.2) but wasn’t successful yet, so I guess I’m missing something. I’m following the documentation and the pendulum example: Automated Index Reduction of DAEs · ModelingToolkit.jl
using ModelingToolkit using OrdinaryDiffEq @parameters t @variables x(t) y(t) a(t) b(t) der = Differential(t) eqs = [der(x) ~ y + a, 0 ~ 3*x - y, der(y) ~ a*x, 0 ~ a/5 - b] @named sys = ODESystem(eqs,t) sys = dae_index_lowering(sys) sys = structural_simplify(sys) u0 = [y => 1.0, x => 3.0, a => 0.0] p =  tspan = (0.0,1.0) prob = ODEProblem(sys, u0, tspan, p) sol = solve(prob, Rodas4())
┌ Warning: Instability detected. Aborting └ @ SciMLBase C:\Users\...\.julia\packages\SciMLBase\n3U0M\src\integrator_interface.jl:351
during solving of the problem and no solution.
I think the problem is in the ODESystem I’m declaring right at the beginning:
@named sys = ODESystem(eqs,t)
In fact the system I want to describe is not an ODE but an DAE.
Also usually I would not need to give start values for all
b. Only one (
y) is a state, so only one of these needs a start value.
I’m not sure if that is the problem or that we (should) have a linear equation system where we need a linear solver to solve several equations at once.