Initialization incomplete for ODESystem on Julia v1.10.10

[oameye]

Describe the bug 🐞

When construct the ODEProblem from a ODESystem, I get the the error:

ERROR: Initialization incomplete. Not all of the state variables of the
DAE system can be determined by the initialization. Missing
variables:

Any[var"⟨σ122⟩"(τ)]

This happens during the construction of the InitializationProblem. Following the runtime with a debugger reveals that everything is fine, i.e., we have the following is empty as it should:

ModelingToolkit.jl/src/problems/initializationproblem.jl

Lines 81 to 83 in f0f412b

	uninit = setdiff(unknowns(sys), unknowns(isys), observables(isys))

until:

ModelingToolkit.jl/src/problems/initializationproblem.jl

Lines 65 to 67 in f0f412b

	if simplify_system
	isys = mtkcompile(isys; fully_determined, split = is_split(sys))
	end

After the mtkcompile one of the observed Equations becomes a parameter_dependency. I.e., before var"⟨σ122⟩"(τ) ~ Initial(var"⟨σ122⟩"(τ)) is in isys.observed, whereas after mtkcompile it is in `parameter dependencies.

Should this happen? Should parameter_dependencies be added to the isempty(uninit) test?

Funnily enough, the system doesn't get "simplified" in Julia 1.10.9. In Julia 1.11.5 the same error occurs

Expected behavior

I expect the the above System compiles to a ODEProblem

Minimal Reproducible Example 👇

Sadly, the system where the errors occurs is fairly large so it is hard to make a smaller reproducible example.

using ModelingToolkit
using QuantumCumulants
using OrdinaryDiffEq

order = 2 #order of the cumulant expansion
@cnumbers N Δ g κ Γ R ν M

hc = FockSpace(:cavity)
ha = NLevelSpace(:atom, 2)

h = hc ⊗ ha

k = Index(h, :k, N, ha)
l = Index(h, :l, N, ha)
m = Index(h, :m, N, ha)
n = Index(h, :n, N, ha)

@qnumbers a::Destroy(h)

σ(i, j, k) = IndexedOperator(Transition(h, :σ, i, j), k)
σn(i, j, k) = NumberedOperator(Transition(h, :σ, i, j), k)

# Define the Hamiltonian
H = -Δ*a'a + g*(Σ(a'*σ(1, 2, k), k) + Σ(a*σ(2, 1, k), k))

J = [a, σ(1, 2, l), σ(2, 1, l), σ(2, 2, l)]
rates = [κ, Γ, R, ν]
ops = [a'*a, σ(2, 2, m)]
eqs = meanfield(ops, H, J; rates = rates, order = order)

# System parameters
N_ = 2e5
Γ_ = 1.0 #Γ=1mHz
Δ_ = 2500Γ_ #Δ=2.5Hz
g_ = 1000Γ_ #g=1Hz
κ_ = 5e6*Γ_ #κ=5kHz
R_ = 1000Γ_ #R=1Hz
ν_ = 1000Γ_ #ν=1Hz

ps = [N, Δ, g, κ, Γ, R, ν]
p0 = [N_, Δ_, g_, κ_, Γ_, R_, ν_]

op1 = σ(1, 2, m)
op2 = σ(2, 1, n)

eqs_c_2 = complete(eqs);

corr3 = CorrelationFunction(op1, op2, eqs_c_2)
corr3_ev = evaluate(corr3, 1, 2; limits = (N=>3));

@named csys3 = ODESystem(corr3_ev)
equations(csys3)[15]

val =ComplexF64[-1.1437146742985789e-7 - 6.3887991023161414e-24im, 0.8635967567168613 + 0.0im, 0.8635967567168613 + 0.0im, 0.8635967567168613 + 0.0im, 3.452543902275645e-7 - 0.00034541449685634095im, 3.4525439022756444e-7 - 0.000345414496856341im, 3.452543902275641e-7 - 0.0003454144968563411im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.00016965965013210335 - 2.727267523449695e-23im, 0.0001696596501321033 - 3.242573084651477e-23im, 0.00016965965013210333 + 4.61495050759439e-24im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.7457993582118818 + 0.0im, 0.7457993582118818 + 3.805486687267425e-24im, 0.7457993582118817 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]
p0_c3 = correlation_p0(corr3_ev, val, ps .=> p0)
u0_c3 = correlation_u0(corr3_ev, val)

dict = merge(Dict(u0_c3), Dict(p0_c3))

prob_c3 = ODEProblem(csys3, dict, (0.0, 0.05))

Error & Stacktrace ⚠️

ERROR: Initialization incomplete. Not all of the state variables of the
DAE system can be determined by the initialization. Missing
variables:

Any[var"⟨σ122⟩"(τ)]

Stacktrace:
  [1] InitializationProblem{…}(sys::System, t::Float64, op::Dict{…}; guesses::Dict{…}, check_length::Bool, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, fully_determined::Bool, check_units::Bool, use_scc::Bool, allow_incomplete::Bool, algebraic_only::Bool, time_dependent_init::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/Documents/ModelingToolkit.jl/src/problems/initializationproblem.jl:91
  [2] InitializationProblem
    @ ~/Documents/ModelingToolkit.jl/src/problems/initializationproblem.jl:20 [inlined]
  [3] #_#1073
    @ ./none:0 [inlined]
  [4] maybe_build_initialization_problem(sys::System, iip::Bool, op::Dict{…}, t::Float64, defs::Dict{…}, guesses::Dict{…}, missing_unknowns::Set{…}; implicit_dae::Bool, time_dependent_init::Bool, u0_constructor::Function, p_constructor::Function, floatT::Type, initialization_eqs::Vector{…}, use_scc::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/Documents/ModelingToolkit.jl/src/systems/problem_utils.jl:1094
  [5] process_SciMLProblem(constructor::Type, sys::System, op::Dict{…}; build_initializeprob::Bool, implicit_dae::Bool, t::Float64, guesses::Dict{…}, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, eval_expression::Bool, eval_module::Module, fully_determined::Bool, check_initialization_units::Bool, u0_eltype::Nothing, tofloat::Bool, u0_constructor::typeof(identity), p_constructor::typeof(identity), check_length::Bool, symbolic_u0::Bool, warn_cyclic_dependency::Bool, circular_dependency_max_cycle_length::Int64, circular_dependency_max_cycles::Int64, substitution_limit::Int64, use_scc::Bool, time_dependent_init::Bool, algebraic_only::Bool, allow_incomplete::Bool, is_initializeprob::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/Documents/ModelingToolkit.jl/src/systems/problem_utils.jl:1352
  [6] (ODEProblem{…})(sys::System, op::Dict{…}, tspan::Tuple{…}; callback::Nothing, check_length::Bool, eval_expression::Bool, expression::Type, eval_module::Module, check_compatibility::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/Documents/ModelingToolkit.jl/src/problems/odeproblem.jl:73
  [7] ODEProblem
    @ ~/Documents/ModelingToolkit.jl/src/problems/odeproblem.jl:65 [inlined]

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

  [35bcea6d] QuantumCumulants v0.4.0
  [961ee093] ModelingToolkit v10.7.0
  [1dea7af3] OrdinaryDiffEq v6.98.0

This enverionment depends on qojulia/QuantumCumulants.jl#243

• Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)

• Output of versioninfo()

Julia Version 1.10.10
Commit 95f30e51f41 (2025-06-27 09:51 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 12 × AMD Ryzen 5 5600X 6-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 10 default, 0 interactive, 5 GC (on 12 virtual cores)
Environment:
  JULIA_EDITOR = code
  JULIA_VSCODE_REPL = 1
  JULIA_NUM_THREADS = 10