Merge pull request #3032 from AayushSabharwal/initial_ode_2

Run initialization on ODEs with superset initial conditions

src/systems/diffeqs/abstractodesystem.jl

@@ -797,6 +797,9 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
     varmap = canonicalize_varmap(varmap)
     varlist = collect(map(unwrap, dvs))
     missingvars = setdiff(varlist, collect(keys(varmap)))
+    setobserved = filter(keys(varmap)) do var
+        has_observed_with_lhs(sys, var) || has_observed_with_lhs(sys, default_toterm(var))
+    end
 
     if eltype(parammap) <: Pair
         parammap = Dict(unwrap(k) => v for (k, v) in todict(parammap))
@@ -810,7 +813,7 @@ 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(setobserved)) &&
          ModelingToolkit.get_tearing_state(sys) !== nothing) ||
         !isempty(initialization_equations(sys))) && t !== nothing
         if eltype(u0map) <: Number

test/components.jl

@@ -47,12 +47,7 @@ sys = structural_simplify(rc_model)
 @test length(equations(sys)) == 1
 check_contract(sys)
 @test !isempty(ModelingToolkit.defaults(sys))
-u0 = [capacitor.v => 0.0
-      capacitor.p.i => 0.0
-      resistor.v => 0.0]
-prob = ODEProblem(sys, u0, (0, 10.0))
-sol = solve(prob, Rodas4())
-check_rc_sol(sol)
+u0 = [capacitor.v => 0.0]
 prob = ODEProblem(sys, u0, (0, 10.0))
 sol = solve(prob, Rodas4())
 check_rc_sol(sol)
@@ -133,9 +128,7 @@ eqs = [connect(source.p, rc_comp.p)
 expand_connections(sys_inner_outer, debug = true)
 sys_inner_outer = structural_simplify(sys_inner_outer)
 @test !isempty(ModelingToolkit.defaults(sys_inner_outer))
-u0 = [rc_comp.capacitor.v => 0.0
-      rc_comp.capacitor.p.i => 0.0
-      rc_comp.resistor.v => 0.0]
+u0 = [rc_comp.capacitor.v => 0.0]
 prob = ODEProblem(sys_inner_outer, u0, (0, 10.0), sparse = true)
 sol_inner_outer = solve(prob, Rodas4())
 @test sol[capacitor.v] ≈ sol_inner_outer[rc_comp.capacitor.v]

test/downstream/inversemodel.jl

@@ -135,7 +135,8 @@ cm = complete(model)
 op = Dict(D(cm.inverse_tank.xT) => 1,
     cm.tank.xc => 0.65)
 tspan = (0.0, 1000.0)
-prob = ODEProblem(ssys, op, tspan)
+# https://github.com/SciML/ModelingToolkit.jl/issues/2786
+prob = ODEProblem(ssys, op, tspan; build_initializeprob = false)
 sol = solve(prob, Rodas5P())
 
 @test SciMLBase.successful_retcode(sol)

test/initial_values.jl

@@ -59,7 +59,9 @@ obs = X2 ~ Γ[1] - X1
 u0 = [X1 => 1.0, X2 => 2.0]
 tspan = (0.0, 1.0)
 ps = [k1 => 1.0, k2 => 5.0]
-@test_nowarn oprob = ODEProblem(osys_m, u0, tspan, ps)
+# Broken since we need both X1 and X2 to initialize Γ but this makes the initialization system
+# overdetermined because parameter initialization isn't in yet
+@test_warn "overdetermined" oprob=ODEProblem(osys_m, u0, tspan, ps)
 
 # Make sure it doesn't error on array variables with unspecified size
 @parameters p::Vector{Real} q[1:3]

test/initializationsystem.jl

@@ -522,3 +522,34 @@ 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.0)
+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[D(Y)][1] == 0.5

test/split_parameters.jl

@@ -82,8 +82,11 @@ eqs = [y ~ src.output.u
 @named sys = ODESystem(eqs, t, vars, []; systems = [int, src])
 s = complete(sys)
 sys = structural_simplify(sys)
-prob = ODEProblem(
+@test_broken ODEProblem(
     sys, [], (0.0, t_end), [s.src.interpolator => Interpolator(x, dt)]; tofloat = false)
+prob = ODEProblem(
+    sys, [], (0.0, t_end), [s.src.interpolator => Interpolator(x, dt)];
+    tofloat = false, build_initializeprob = false)
 sol = solve(prob, ImplicitEuler());
 @test sol.retcode == ReturnCode.Success
 @test sol[y][end] == x[end]