`solve` fails with `InitialFailure` : overdetermined system

[harisorgn]

Describe the bug 🐞

solve fails with InitialFailure because it treats my system as overdetermined. System is technically not overdetermined, it just has an state (jcn(t)) that only appears in an algebraic equation and it is successfully substituted/eliminated during structural_simplify.

Expected behavior

A clear and concise description of what you expected to happen.

Minimal Reproducible Example 👇

IIUC generate_initializesystem treats the default value of jcn(t)=0 as one equation in the InitialSystem and its algebraic equation jcn ~ -hx as a second equation, both involving a single state, thus overdetermined system? If I remove the default value of jcn(t) it works fine.

Possibly related to #3033

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

@variables x(t)=1 jcn(t)=0 [input = true]
@parameters k=1 h=2

eqs = [
    D(x) ~ -k*x + jcn,
    jcn ~ -h*x
]

@named sys = System(eqs, t)
ss = structural_simplify(sys)
prob = ODEProblem(ss, [], (0, 10)) # errors with ;fully_determined=true
sol = solve(prob, Tsit5())

Error & Stacktrace ⚠️

 sol = solve(prob, Tsit5()) # retcode: InitialFailure

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

[961ee093] ModelingToolkit v9.39.1
⌃ [1dea7af3] OrdinaryDiffEq v6.87.0

• Output of versioninfo()

Julia Version 1.10.0
Commit 3120989f39 (2023-12-25 18:01 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 20 × 13th Gen Intel(R) Core(TM) i5-13600KF
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, goldmont)
  Threads: 1 on 20 virtual cores
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS =

[hersle]

Isn't the initial value of jac uniquely determined by the 2nd equation when the initial value of x is specified? Then I would expect this behavior when a default value for jac is provided, too.

[harisorgn]

Isn't the initial value of jac uniquely determined by the 2nd equation when the initial value of x is specified?

Yep, this would be the expected behaviour for me too. For context, this is part of a higher level interface where jcn-type variables can have input equations like the jcn ~ -h*x above, but such equations depend on other Systems that the current System connects to. If there is no such connection made, then the fallback is to use the default value jcn(t)=0.

[ChrisRackauckas]

Yes, there is nothing wrong here. It fails initialization because it's overdetermined and 0 is not the solution, but you're saying 0 has to be the solution. The working code is of course:

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

@variables x(t)=1 jcn(t) [input = true]
@parameters k=1 h=2

eqs = [
    D(x) ~ -k*x + jcn,
    jcn ~ -h*x
]

@named sys = System(eqs, t)
ss = structural_simplify(sys)
prob = ODEProblem(ss, [], (0, 10)) # errors with ;fully_determined=true
sol = solve(prob, Tsit5())

Do not mix guesses with chosen defaults.

[ChrisRackauckas]

If there is no such connection made, then the fallback is to use the default value jcn(t)=0.

That's not how the default works. It is always a default equation, jcn ~ 0 at initialization. I recommend just connecting a constant term to close the system, that's required anyways to end up with a system that fully connects all inputs (all inputs being connected is a standard correctness check)

[harisorgn]

The working code is of course

Yep, this was my fix. I just opened the issue for others, because I was not certain if this is expected behaviour or bug.

... but you're saying 0 has to be the solution ... Do not mix guesses with chosen defaults.

Got it, I was treating it as a guess.

[ChrisRackauckas]

No, a guess is different. jcn(t) [input = true, guess=0]. If you give a value/equation, MTK will make sure it's always satisfied. If it's supposed to be a guess, then it should be explicitly marked as a guess. This is one of the v9 changes.