DAE simulation fails even when all derivatives are zero

Specific Domains Modelling & Simulations
1

Hello!

I have this system:

using ModelingToolkit
using ModelingToolkit: t_nounits, D_nounits
using OrdinaryDiffEqRosenbrock, OrdinaryDiffEqNonlinearSolve

@variables begin
    cation1m(t_nounits), [guess=1.0]
    anion1m(t_nounits), [guess=1.0]
    salt1M(t_nounits), [guess=1.0]
    anion2m(t_nounits), [guess=1.0]
    salt2M(t_nounits), [guess=1.0]
    totalcat(t_nounits), [guess=2.0]
    totalan1(t_nounits), [guess=2.0]
    totalan2(t_nounits), [guess=2.0]
    extra_an1(t_nounits), [guess=0.0]
    extra_an2(t_nounits), [guess=0.0]
end
@parameters begin
    maw1=2.0
    maw2=3.0
end
eqs = [
    W ~ salt1M / maw1 + salt2M / maw2
    cation1m * W ~ salt1M + salt2M
    anion1m * W ~ salt1M
    anion2m * W ~ salt2M
    totalcat ~ salt1M + salt2M
    totalan1 ~ salt1M + extra_an1
    totalan2 ~ salt2M + extra_an2
    0 ~ cation1m - anion1m - anion2m
    D_nounits(totalcat) ~ 0.0
    D_nounits(totalan1) ~ 0.0
    D_nounits(totalan2) ~ 0.0
]
@mtkcompile aqt = System(eqs, t_nounits)

prob = ODEProblem(aqt,
    [
        aqt.totalcat => 2.0,
        aqt.totalan1 => 1.5,
        aqt.totalan2 => 2.0,
    ],
    (0.0, 1.0))

sol = solve(prob, Rosenbrock23()) # at t=0, dt was forced below floating point epsilon...

The code above gives the following error:

┌ Warning: At t=0.0, dt was forced below floating point epsilon 5.0e-324, and step error estimate = NaN. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase ~/.julia/packages/SciMLBase/YE7xF/src/integrator_interface.jl:657
retcode: Unstable
Interpolation: specialized 2nd order "free" stiffness-aware interpolation
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Vector{Float64}}:
 [0.8433629863955399, 1.2570782481872695, 1.1143800700692865, 2.0, 1.5, 2.0]

My question is, since the system was able to initialize correctly, and the derivatives are all zero, why wasn’t the simulation able complete? Thanks for any suggestions you can offer!

2

Rosenbrock23 is not fully stable for DAEs. Does Rodas5P work?

3

Thanks for your response! Rodas5P doesn’t work either, however.

sol = solve(prob, Rodas5P())

┌ Warning: At t=0.0, dt was forced below floating point epsilon 5.0e-324, and step error estimate = NaN. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase ~/.julia/packages/SciMLBase/YE7xF/src/integrator_interface.jl:657
retcode: Unstable
Interpolation: specialized 4th (Rodas6P = 5th) order "free" stiffness-aware interpolation
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Vector{Float64}}:
 [0.8433629863955399, 1.2570782481872695, 1.1143800700692865, 2.0, 1.5, 2.0]
4

When letting the solver decide the solver by itself it works for me:

using DifferentialEquations

...

sol = solve(prob)

I’m not sure why the problem as defined is so difficult to solve. I haven’t checked the index of the DAE but maybe that’s the problem? Also, did you set the derivatives to zero to find the steady-state solution or initial condition? Maybe we can discuss some different methods to do that.

1 Like
5

That works, thank you!

The problem above is just a MWP version of a larger system found here. Ultimately, the goal isn’t to find the steady state, but to use it as a subsystem that connects into dynamics from another system. Unfortunately it just gets harder to solve as the system gets larger, so I’d be very interested in suggestions regarding why this small system is so hard to solve, or a different way to formulate it that would be easier. Thanks for your thoughts!

6

Here’s a version of the problem that also fails with the default polyalgorithm:

using ModelingToolkit
using ModelingToolkit: t_nounits, D_nounits
using OrdinaryDiffEq

@variables begin
    cation1m(t_nounits), [guess=1.0]
    anion1m(t_nounits), [guess=1.0]
    salt1M(t_nounits), [guess=1.0]
    anion2m(t_nounits), [guess=1.0]
    salt2M(t_nounits), [guess=1.0]
    W(t_nounits), [guess=1.0]
    totalcat(t_nounits), [guess=2.0]
    totalan1(t_nounits), [guess=2.0]
    totalan2(t_nounits), [guess=2.0]
    solid1(t_nounits), [guess=0.0]
    solid2(t_nounits), [guess=0.0]
end
@parameters begin
    maw1=2.0
    maw2=3.0
end
eqs = [
    W ~ salt1M / maw1 + salt2M / maw2
    cation1m * W ~ salt1M + salt2M
    anion1m * W ~ salt1M
    anion2m * W ~ salt2M
    totalcat ~ salt1M + salt2M
    totalan1 ~ salt1M + solid1 + 2solid2
    totalan2 ~ salt2M + 2solid1 + solid2
    0 ~ cation1m - anion1m - anion2m
    D_nounits(totalcat) ~ 0.0
    D_nounits(totalan1) ~ 0.0
    D_nounits(totalan2) ~ 0.0
]
@mtkcompile aqt = System(eqs, t_nounits)

prob = ODEProblem(aqt,
    [
        aqt.totalcat => 2.0,
        aqt.totalan1 => 1.5,
        aqt.totalan2 => 2.0,
    ],
    (0.0, 1.0))

sol = solve(prob) # ┌ Warning: At t=0.0, dt was forced below floating point epsilon 5.0e-324, and step error estimate = NaN.
7

Looking at your system of equations, I can see that the 2 equations with cation1m are actually the same. You need to swap out one of these equations with some new information for the numerical solver to work.

I wrote a list of debugging tips for MTK here: Debugging difficult stiff ODE/DAE models · ModelingToolkit Course It’s a bit out of date since the new releases of MTK, but the ideas are still valid.

Hopefully something there may be of some help. This stuff is always very challenging, hope you can find a small system that works and helps you move the ball forward.

3 Likes
8

Thanks, I didn’t realize that! I guess I had kind of assumed (naively?) that the simplification process would have caught that case and removed one of the equations. But now I know that it doesn’t!