Supplying `u0` to `ODEProblem` does not work in some cases

[bradcarman]

If I use InitializationProblem to pre-calculate a u0 for a simple problem, this works without incident...

@mtkmodel Simple begin
    @variables begin
        x(t) = 0
    end
    @equations begin
        D(x) ~ 1
    end
end

@mtkbuild sys = Simple()
initprob = ModelingToolkit.InitializationProblem(sys, 0.0)
initsol = solve(initprob)
u0 = initsol[unknowns(sys)]
prob = ODEProblem(sys, u0, (0, 0))  #OK, no warning

However, for a more complex problem, for some reason I do not get the same behavior...

function System()
    pars = @parameters begin
        rho_0 = 1000
        beta = 2e9
        A = 0.1
        m = 100
        L = 1
        p_s = 100e5
        p_r = 10e5
        C = 2.7*0.5
        c = 1000
        A_p = 0.00094
    end   
    vars = @variables begin
        p_1(t) = p_s
        p_2(t) = p_r
        x(t)=0
        dx(t)=0
        ddx(t), [guess=0]
        rho_1(t), [guess=rho_0]
        rho_2(t), [guess=rho_0]
        drho_1(t), [guess=0]
        drho_2(t), [guess=0]
        dm_1(t), [guess=0]
        dm_2(t), [guess=0]
    end
    
    # let -----
    u_1 = dm_1/(rho_0*A_p)
    u_2 = dm_2/(rho_0*A_p)

    eqs = [
        D(x) ~ dx
        D(dx) ~ ddx
        D(rho_1) ~ drho_1
        D(rho_2) ~ drho_2
        +dm_1 ~ drho_1*(L+x)*A + rho_1*dx*A
        -dm_2 ~ drho_2*(L-x)*A - rho_2*dx*A
        rho_1 ~ rho_0*(1 + p_1/beta)
        rho_2 ~ rho_0*(1 + p_2/beta)
        m*ddx ~ (p_1 - p_2)*A - c*dx
        (p_s - p_1) ~ C*rho_0*(u_1)*abs(u_1)
        (p_2 - p_r) ~ C*rho_0*(u_2)*abs(u_2)
    ]

    return ODESystem(eqs, t, vars, pars; name)
end

@mtkbuild sys = System()
initprob = ModelingToolkit.InitializationProblem(sys, 0.0)
initsol = solve(initprob)
u0 = initsol[unknowns(sys)]
prob = ODEProblem(sys, u0, (0, 0)) # Warning: Initialization system is overdetermined. 4 equations for 0 unknowns.

Note: my system is not overdetermined! I can see this because for one, the InitializationProblem build and solves just fine. Also I can simply remove u0 without incident...

julia> prob = ODEProblem(sys, [], (0, 0))
ODEProblem with uType Vector{Float64} and tType Int64. In-place: true
timespan: (0, 0)
u0: 6-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

[ChrisRackauckas]

Yeah this looks like an interface bug.

[AayushSabharwal]

This isn't an interface bug. In the first example, there is a single variable x, and the initial equation x ~ 0 from the default. Passing this solution to ODEProblem doesn't emit a warning because there are no algebraic variables, and all differential variables have an initial condition.

In your second example, there are differential and algebraic variables. ModelingToolkit.InitializationProblem(sys, 0.0) is fully determined because:

• Two of the differential variables (x and dx) have defaults which act as initial equations
• The other two differential variables (rho_1 and rho_2) can be solved from the observed equations for p_1 and p_2, where p_1 and p_2 have defaults
• The remaining four variables in the simplified system (drho_1, drho_2, dm_1, dm_2) are algebraic and solvable from the corresponding equations.

When providing the solution from initprob to ODEProblem, the ODEProblem thinks it has 8 initial conditions, four of which are for algebraic variables. It has to ensure that the algebraic equations are satisfied, so it makes an InitializationProblem from the 8 initial conditions + 4 algebraic equations and thus is overdetermined. Since the 8 initial conditions are consistent, it will successfully solve. However, it is still an overdetermined system.

[AayushSabharwal]

As a contrived example, consider the system

D(x) ~ x
x^3 + y^3 ~ 3

Where x => 1.0 is a default. The initialization system is:

x ~ 1.0
x^3 + y^3 ~ 3

Which solves to [x => 1.0, y => cbrt(2)]. If you provide this initial condition to ODEProblem, it will still create an InitializationProblem to verify the algebraic equations. Thus, it creates the system

x ~ 1.0
y ~ cbrt(2)
x^3 + y^3 ~ 3

which is a valid initialization, but overdetermined nonetheless.

[bradcarman]

My understanding of how ModelingToolkit works, is that if you supply u0 as a vector of numbers, then this is treated as the provided u0 for the ODEProblem and the initialization system is bypassed. It seems in my first example, this is the behavior I get. Furthermore, in the example you provide, I also get this behavior...

@mtkmodel Simple begin
    @variables begin
        x(t) = 1
        y(t), [guess=cbrt(2)]
    end
    @equations begin
        D(x) ~ x
        x^3 + y^3 ~ 3
    end
end

@mtkbuild sys = Simple()
initprob = ModelingToolkit.InitializationProblem(sys, 0.0)
initsol = solve(initprob)
u0 = initsol[unknowns(sys)]
prob = ODEProblem(sys, u0, (0, 0))  #OK, no warning

So I'm confused, why in my 2nd example, do I not get this same behavior? I'm supplying a Vector{Float64} as u0, there should be no Initialization System invoked, just like in the example above.

[AayushSabharwal]

My understanding of how ModelingToolkit works, is that if you supply u0 as a vector of numbers, then this is treated as the provided u0 for the ODEProblem and the initialization system is bypassed.

This isn't always the case. A vector u0 is treated as unknowns(sys) .=> u0. This alone is sufficient to fully determine initial values for all unknowns, and hence the initialization system sometimes doesn't need to be built. In the @mtkmodel Simple example, this is what happens. You gave 2 initial conditions for 2 variables, and there were no additional defaults in the mix (the default for x is overridden by the initial condition for x).

In your second example, in addition to the 8 initial conditions from u0 there are additional conditions from defaults of observed variables. Namely: the defaults for p_1 and p_2. To ensure these are respected, MTK creates the initialization problem which ends up being overdetermined. MTK basically says "I know p_1(t) ~ (beta*(-rho_0 + rho_2(t))) / rho_0 from the observed equations, and the user has given a default of p_1(t) => p_s. I need to make sure both of these are satisfied given the rest of the initial conditions". To summarize, there are 8 initial conditions from u0, 3 from observed(sys), 2 from the defaults for p_1 and p_2 and 4 from the algebraic equations totalling to 17 conditions. If you look at prob.f.initialization_data.f.sys, it has 6 equations and 11 observed equations - the same 17 conditions.

I apologize for the confusion from my example, I was a little mistaken and should have run some code rather than trying to emulate what MTK would do from memory. However, I do feel we should have that behavior, otherwise the user thinks the initial condition they provided is respected and BrownBasicInit will subsequently change it when actually solving.

@ChrisRackauckas In this example:

julia> @mtkmodel Simple begin
           @variables begin
               x(t) = 1, [guess = 1]
               y(t) = cbrt(2), [guess=cbrt(2)]
           end
           @equations begin
               D(x) ~ x
               x^3 + y^3 ~ 3
           end
       end
julia> @mtkbuild sys = Simple()
julia> prob = ODEProblem(sys, [sys.y => cbrt(3)], (0.0, 1.0))  #OK, no warning
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 2-element Vector{Float64}:
 1.0
 1.4422495703074083

MTK is more than happy to create the system without an initialization problem. Creating the integrator:

julia> integ = init(prob)
t: 0.0
u: 2-element Vector{Float64}:
 1.0
 1.2599210499101785

BrownBasic corrects the algebraic variable. However, it is entirely possible that the initial condition for x is wrong and should be 0.0 instead.

[AayushSabharwal]

If you want to unconditionally bypass the initialization problem, you need to pass build_initializeprob = false to ODEProblem.

[AayushSabharwal]

Here's the condition for creating an initialization problem, paraphrased from the if condition in the code

(
    (
        ("we are constructing a DAEProblem" OR "there are initial conditions for observed variables" OR "there are unknowns that don't have an initial condition" OR "there are parameters to solve for" OR "the initial condition for an unknown is a symbolic expression")
        AND
        "the system was structurally simplified"
    )
    OR
    "the system has initialization_eqs"
)
AND
"the system is time-dependent and the user has provided an initial time to the problem constructor"

[ChrisRackauckas]

BrownBasic corrects the algebraic variable. However, it is entirely possible that the initial condition for x is wrong and should be 0.0 instead.

That case shouldn't bypass initialization. MTK shouldn't do that for any case with algebraic variables?

[AayushSabharwal]

Yeah, I was confused to find it bypassed initialization. Will fix ASAP.

[bradcarman]

OK, thanks for explaining the build_initializeprob = false keyword. Note: this is not documented.

So, now I'm even more confused. In the code below, I provide initial equations init_eqs that are not possible and are overdetermined for this system. When I use InitializationProblem I get the expected warning, and when I attempt to solve, I get a Failure retcode. This is all expected. However, when I supply these initial equations to the ODEProblem I get no warning. So now I'm purposely defining an overdetermined system, and I don't get a warning. This is very confusing. I think ODEProblem needs to be more transparent about what's happening with initialization. In this case the user is expecting to get [sys.x => 1, sys.y => cbrt(3)] but they get something different.

@mtkmodel Simple begin
    @variables begin
        x(t), [guess=1]
        y(t), [guess=1]
    end
    @equations begin
        D(x) ~ x
        x^3 + y^3 ~ 3
    end
end

@mtkbuild sys = Simple()

init_eqs = [sys.x => 1, sys.y => cbrt(3)]

initprob = ModelingToolkit.InitializationProblem(sys, 0.0, init_eqs) # Warning: Initialization system is overdetermined
initsol = solve(initprob) # retcode: Failure

prob = ODEProblem(sys, init_eqs, (0, 0.1))  #OK, no warning

[AayushSabharwal]

Note: this is not documented.

Would adding this and other such keyword arguments to the ODEProblem docstring help?

In this case the user is expecting to get [sys.x => 1, sys.y => cbrt(3)] but they get something different.

Yes, this is the behavior Chris and I were discussing should be fixed. Currently, MTK thinks "I have all the initial conditions, don't need an initialization problem" but that isn't the case and it should build an initprob because the system has algebraic equations. It's also the behavior I expected in the example I gave here. After the fix, providing that initial condition will InitialFailure when the ODEProblem is solved.

[bradcarman]

Yes, adding the keyword args that ODEProblem accept to the doc string would be really helpful. Thanks!

[ChrisRackauckas]

Would adding this and other such keyword arguments to the ODEProblem docstring help?

We shouldn't really document that. It's an internal that would lead to correctness issues, it shouldn't be exposed to users since it would get easily abused.