Determined complete system raises an overdetermined warn when creating ODEProblem

[pmc4]

Describe the bug 🐞

I have created a system with structural_simplify that has 2 equations, 2 unknowns and 8 observed variables. The system is marked as complete and it's a determined system of equations. When I create the problem with ODEProblem, a warning stating that the initialization system is overdetermined appears, saying that there are 7 equations for 0 unknowns.

I can still solve the problem and get a success retcode, though.

Expected behavior

I should be able to pass a complete determined system to the ODEProblem constructor and create the problem without warnings.

Minimal Reproducible Example 👇

I have simplified a lot the structure of my models, but in essence is what I've written below. I know I'm using @independent_variables x instead of t_nounits.

using ModelingToolkit, DifferentialEquations, Integrals

@independent_variables x
D = Differential(x)

rhoA(z) = z^4
rhoB(z) = z^3
PB(z) = z^4

# Define models and solve
@mtkmodel ModelA begin
    @variables begin
        z(x)
        ρ(x)
        P(x)
        ρ̇(x)
    end
    @equations begin
        ρ ~ rhoA(z)
        P ~ ρ / 3
        ρ̇ ~ 0
    end
end


@mtkmodel ModelB begin
    @variables begin
        z(x)
        ρ(x)
        P(x)
        ρ̇(x)
    end
    @equations begin
        ρ ~ rhoB(z)
        P ~ PB(z)
        ρ̇ ~ (ρ - 3*P) / x
    end
end


function Model()
    z0 = 1.0

    # Model A
    ρ0A = rhoA(z0)
    P0A = ρ0A / 3
    z0A = z0
    @named A = ModelA(ρ=ρ0A, P=P0A, z=z0A)

    # Model B
    ρ0B = rhoB(z0)
    P0B = PB(z0)
    z0B = z0
    @named B = ModelB(ρ=ρ0B, P=P0B, z=z0B)

    system = [A, B]
    ρtot0 = ρ0A + ρ0B

    @variables z(x) = z0 ρtot(x)

    eqs = [
        ρtot ~ A.ρ + B.ρ
        D(ρtot) ~ A.ρ̇ + B.ρ̇
        A.z ~ z
        B.z ~ z
    ]

    @named _model = ODESystem(eqs, x)
    @named model = compose(_model, system)

    return model, ρtot0
end


model, ρtot0 = Model()
sys = structural_simplify(model)


xspan = (1e-2, 1e2)
prob = ODEProblem(sys, [sys.ρtot => ρtot0], xspan, []; warn_initialize_determined = true)
sol = solve(prob)

Error & Stacktrace ⚠️

┌ Warning: Initialization system is overdetermined. 7 equations for 0 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/41wGH/src/systems/diffeqs/abstractodesystem.jl:1358

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

Status `~/software/test/Project.toml`
  [0c46a032] DifferentialEquations v7.15.0
  [de52edbc] Integrals v4.5.0
  [961ee093] ModelingToolkit v9.61.0

• Output of versioninfo()

Julia Version 1.11.2
Commit 5e9a32e7af2 (2024-12-01 20:02 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 20 × 12th Gen Intel(R) Core(TM) i7-12700H
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, alderlake)
Threads: 1 default, 0 interactive, 1 GC (on 20 virtual cores)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 

[AayushSabharwal]

This is an overdetermined system.

You have one algebraic equation:

0 ~ -ρtot(x) + B₊ρ(x) + A₊ρ(x)

Seven defaults:

Dict{Any, Any} with 7 entries:
  z(x)   => 1.0
  A₊ρ(x) => 1.0
  B₊P(x) => 1.0
  A₊P(x) => 0.333333
  B₊ρ(x) => 1.0
  A₊z(x) => 1.0
  B₊z(x) => 1.0

Eight observed equations:

8-element Vector{Equation}:
 B₊P(x) ~ B₊z(x)^4
 B₊ρ(x) ~ B₊z(x)^3
 A₊z(x) ~ B₊z(x)
 A₊ρ̇(x) ~ 0.0
 B₊ρ̇(x) ~ (B₊ρ(x) - 3B₊P(x)) / x
 A₊ρ(x) ~ A₊z(x)^4
 z(x) ~ A₊z(x)
 A₊P(x) ~ (1//3)*A₊ρ(x)

And one initial condition

[sys.ρtot => ρtot0]

Which totals to 17 initial conditions. Taking a look at the initialization system:

julia> isys = prob.f.initialization_data.initializeprob.f.sys
Model model:
Equations (7):
  7 standard: see equations(model)
Parameters (1): see parameters(model)
  x
Observed (10): see observed(model)

7 equations + 10 observed == 17 total

Note:

julia> full_equations(isys)
ERROR: tearing failed because the system is singular

The system is singular because a lot of your conditions are redundant. For example,

ρtot(x) ~ 2.0 # from the initial condition provided to `ODEProblem`
A₊ρ(x) ~ 1.0 # from defaults
B₊P(x) ~ 1.0 # from defaults
0 ~ -ρtot(x) + B₊ρ(x) + A₊ρ(x) # algebraic equation

4 equations, 3 variables. Given any 3 of these equations, the fourth one is automatically true. The reason the system can still be solved is that the overdetermined conditions are still consistent. For example, if you had given B.P a default of 2.0, the solve would fail with InitialFailure.

[pmc4]

You're completely right. In fact my written equations can be transformed into dz/dx and the only initial condition I need is z0 for a given x0. I got confused and since I wasn't getting that warning at first, I thought it had something to do with an update.

Thanks a lot for figuring it out and I'm sorry for the bothering!