Run initialization on ODEs with superset initial conditions

InitializationProblem is not always built, however it missed a case. If you have an ODE and you define all of the initial conditions, then it does not need an initialization problem. But, if you simplify down to an ODE and give initial conditions for all of the original variables, including some simplified out, then you have defined initial conditions on some observed quantities and that needs to be accounted for by an initialization process.

This fixes #2619

Author: ChrisRackauckas

src/systems/diffeqs/abstractodesystem.jl

@@ -797,6 +797,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
     varmap = canonicalize_varmap(varmap)
     varlist = collect(map(unwrap, dvs))
     missingvars = setdiff(varlist, collect(keys(varmap)))
+    setboserved = setdiff(collect(keys(varmap)), varlist)
 
     if eltype(parammap) <: Pair
         parammap = Dict(unwrap(k) => v for (k, v) in todict(parammap))
@@ -810,7 +811,8 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
 
     # ModelingToolkit.get_tearing_state(sys) !== nothing => Requires structural_simplify first
     if sys isa ODESystem && build_initializeprob &&
-       (((implicit_dae || !isempty(missingvars)) &&
+       (((implicit_dae || !isempty(missingvars) || !isempty(setboserved)) &&
+         all(isequal(Continuous()), ci.var_domain) &&
          ModelingToolkit.get_tearing_state(sys) !== nothing) ||
         !isempty(initialization_equations(sys))) && t !== nothing
         if eltype(u0map) <: Number

test/initializationsystem.jl

@@ -522,3 +522,32 @@ end
           structural_simplify
     @test_nowarn ODEProblem(sys, [D(x) => 1.0], (0.0, 1.0), [])
 end
+
+# https://github.com/SciML/ModelingToolkit.jl/issues/2619
+@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)
+
+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]
+
+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)
+
+@test solve(oprob_1st_order_1, Rosenbrock23()).retcode == SciMLBase.ReturnCode.InitialFailure
+@test solve(oprob_1st_order_2, Rosenbrock23())[Y][1] == 2.0
+@test solve(oprob_2nd_order_1, Rosenbrock23()).retcode == SciMLBase.ReturnCode.InitialFailure
+sol = solve(oprob_2nd_order_2, Rosenbrock23()) # retcode: Success
+@test sol[Y][1] == 2.0
+@test sol[1][2] == 0.5