Case of unstable initial condition finding with BrownFullBasicInit

[AnHeuermann]

This issue is from a discussion on discourse.julialang.org about ModelingToolkit.jl and that I used wrong start values for my ODE problem.

In short

I'm using the wrong start values u0 for x and y and during simulation I'll get an error:

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

At that time I could not find the source of the error.

---

My suggestion would be to improve the error message so that a user gets an error like

Error: The initialization problem is inconsistent due to the following equation: 0 != 8 = 3*x - y
Aborting

(Error message in parts taken from OpenModelica).

I would argue for an error over a warning because then the user needs to fix the model and can't simply ignore the problem only the see it pop up in the next tool again.

[AnHeuermann]

CC @ChrisRackauckas

[ChrisRackauckas]

That's not the issue. As mentioned on Discourse, this solver can automatically find initial conditions. The issue is that if you check the return:

(solve(prob, Rodas4()))[1] = [3.0, 1.0, -Inf]

for some reason the initial condition finder failed and return a -Inf. So it's an error there.

Did you try the Shampine initialization?

[ChrisRackauckas]

We should add a warn and exit if the initialization fails though.

[AnHeuermann]

Did you try the Shampine initialization?

Not yet, but when I try it now it simulates, but the solution seems to be wrong and I don't seem to have the control on which start values are used.

u0 = [y => 1.0,
      x => 3.0,
      a => 0.0]
prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, Rodas4(),  initializealg=ShampineCollocationInit())

I'll get sol(0) = [2.949931674070046, 0.8497950222101376, -50.91812095216428]=[x,y,a]. This violates equation y = 3*x.

The correct solution for this y=0.85 should be x=0.28, a=-0.94 (According to OpenModelica).

[AnHeuermann]

The wrong solution is probably because of sys = dae_index_lowering(sys) and not an issue of OrdinaryDiffEq.jl but SciML/ModelingToolkit.jl. It is simply that an important equation y=3*x from the original system is missing in the reduced system.

The dummy-derivate method would be one way to solve this. I can open a ticket about this on SciML/ModelingToolkit.jl.

---

EDIT: Ah I see @kabdelhak already commented this on https://github.com/SciML/ModelingToolkit.jl/issues/988#issuecomment-935762807 yesterday.

[ChrisRackauckas]

Yes, that looks like it could be what's going on. Nice snooping.