Inconsistency for how to give intiial conditions to DAEs with higher order derivatives.

[TorkelE]

Normally, one can/should give initial values/guesses for variables occurring only in the algebraic equations parts of a DAE. However, when one does this for a DAE with higher-order derivatives there are weird behaviours:

using OrdinaryDiffEq, ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

# Creates systems.
@parameters k1 k2 ω
@variables X(t) Y(t)
eqs_1st_order = [D(Y) + Y - ω ~ 0,
       X + k1 ~ Y + k2]
eqs_2nd_order = [D(D(Y)) + 2ω*D(Y) + (ω^2)*Y ~ 0,
        X + k1 ~ Y + k2]
@mtkbuild sys_1st_order = ODESystem(eqs_1st_order, t)
@mtkbuild sys_2nd_order = ODESystem(eqs_2nd_order, t)

# Creates initial conditions and parameters.
u0_1st_order_1 = [X => 1.0, Y => 2.0]
u0_1st_order_2 = [Y => 2.0]
u0_2nd_order_1 = [X => 1.0, Y => 2.0, D(Y) => 0.5]
u0_2nd_order_2 = [Y => 2.0, D(Y) => 0.5]
tspan = (0.0, 10.)
ps = [ω => 0.5, k1 => 2.0, k2 => 3.0]

# Creates ODEProblems.
oprob_1st_order_1 = ODEProblem(sys_1st_order, u0_1st_order_1, tspan, ps)
oprob_1st_order_2 = ODEProblem(sys_1st_order, u0_1st_order_2, tspan, ps)
oprob_2nd_order_1 = ODEProblem(sys_2nd_order, u0_2nd_order_1, tspan, ps) # gives sys_2nd_order
oprob_2nd_order_2 = ODEProblem(sys_2nd_order, u0_2nd_order_2, tspan, ps)

# Solves ODEProblems.
solve(oprob_1st_order_1, Rosenbrock23()) # retcode: Success
solve(oprob_1st_order_2, Rosenbrock23()) # retcode: Success
solve(oprob_2nd_order_1, Rosenbrock23()) # retcode: InitialFailure
solve(oprob_2nd_order_2, Rosenbrock23()) # retcode: Success

[ChrisRackauckas]

What is the weird behavior?

[TorkelE]

Here, if you provide an initial guess for X (the variable that only appears in the algebraic equation) when attempting to solve the 2nd-order DAE, you get an error. However, it is possible to give such a guess to the DAE with first order differentials only.

(This might make sense for someone more involved, but encountered it when running various tests, and figured that it seemed weird)

[ChrisRackauckas]

That's not a guess, it's a condition, and the 3 conditions together are not satisfiable. If you want to supply a guess, you have to make it a guess. See guesses in https://docs.sciml.ai/ModelingToolkit/stable/tutorials/initialization/

[TorkelE]

Right, that makes sense. But why do I then not get the same problem for the first case? I.e.

u0_1st_order_1 = [X => 1.0, Y => 2.0]
oprob_1st_order_1 = ODEProblem(sys_1st_order, u0_1st_order_1, tspan, ps)
sol = solve(oprob_1st_order_1, Rosenbrock23()) # retcode: Success

then when I check

sol[X][1]

I get 3.0 (which is different from the initial condition I gave)

[ChrisRackauckas]

That's a bug because it simplifies down to a pure ODE and then doesn't run the initialization consistency check that DAEs do. Maybe I can address that in https://github.com/SciML/ModelingToolkit.jl/pull/2591