Fix guess and initial condition length checking

Author: ChrisRackauckas

src/systems/abstractsystem.jl

@@ -2268,14 +2268,14 @@ function UnPack.unpack(sys::ModelingToolkit.AbstractSystem, ::Val{p}) where {p}
 end
 
 """
-    missing_variable_defaults(sys::AbstractSystem, default = 0.0)
+    missing_variable_defaults(sys::AbstractSystem, default = 0.0; subset = unknowns(sys))
 
 returns a `Vector{Pair}` of variables set to `default` which are missing from `get_defaults(sys)`.  The `default` argument can be a single value or vector to set the missing defaults respectively.
 """
-function missing_variable_defaults(sys::AbstractSystem, default = 0.0)
+function missing_variable_defaults(sys::AbstractSystem, default = 0.0; subset = unknowns(sys))
     varmap = get_defaults(sys)
     varmap = Dict(Symbolics.diff2term(value(k)) => value(varmap[k]) for k in keys(varmap))
-    missingvars = setdiff(unknowns(sys), keys(varmap))
+    missingvars = setdiff(subset, keys(varmap))
     ds = Pair[]
 
     n = length(missingvars)

src/systems/diffeqs/abstractodesystem.jl

@@ -848,18 +848,34 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
         tofloat = true,
         symbolic_u0 = false,
         u0_constructor = identity,
+        guesses = Dict(),
+        warn_initialize_determined = true,
         kwargs...)
     eqs = equations(sys)
     dvs = unknowns(sys)
     ps = full_parameters(sys)
     iv = get_iv(sys)
 
+    initializeprob =  ModelingToolkit.InitializationProblem(sys, u0map, parammap; guesses, warn_initialize_determined)
+    initializeprobmap = getu(initializeprob, unknowns(sys))
+
+    # Append zeros to the variables which are determined by the initialization system
+    # This essentially bypasses the check for if initial conditions are defined for DAEs
+    # since they will be checked in the initialization problem's construction
+    # TODO: make check for if a DAE cheaper than calculating the mass matrix a second time!
+    if implicit_dae || calculate_massmatrix(sys) !== I
+        zerovars = setdiff(unknowns(sys),defaults(sys)) .=> 0.0
+        trueinit =  identity.([zerovars;u0map])
+    else
+        trueinit = u0map
+    end
+
     if has_index_cache(sys) && get_index_cache(sys) !== nothing
-        u0, defs = get_u0(sys, u0map, parammap; symbolic_u0)
+        u0, defs = get_u0(sys, trueinit, parammap; symbolic_u0)
         p = MTKParameters(sys, parammap)
     else
         u0, p, defs = get_u0_p(sys,
-            u0map,
+            trueinit,
             parammap;
             tofloat,
             use_union,
@@ -886,9 +902,6 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
 
     check_eqs_u0(eqs, dvs, u0; kwargs...)
 
-    initializeprob =  ModelingToolkit.InitializationProblem(sys, u0map, parammap)
-    initializeprobmap = getu(initializeprob, unknowns(sys))
-
     f = constructor(sys, dvs, ps, u0; ddvs = ddvs, tgrad = tgrad, jac = jac,
         checkbounds = checkbounds, p = p,
         linenumbers = linenumbers, parallel = parallel, simplify = simplify,
@@ -998,13 +1011,14 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map =
         parammap = DiffEqBase.NullParameters();
         callback = nothing,
         check_length = true,
+        warn_initialize_determined = true,
         kwargs...) where {iip, specialize}
     if !iscomplete(sys)
         error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`")
     end
     f, u0, p = process_DEProblem(ODEFunction{iip, specialize}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
-        check_length, kwargs...)
+        check_length, warn_initialize_determined, kwargs...)
     cbs = process_events(sys; callback, kwargs...)
     inits = []
     if has_discrete_subsystems(sys) && (dss = get_discrete_subsystems(sys)) !== nothing
@@ -1069,13 +1083,14 @@ end
 
 function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan,
         parammap = DiffEqBase.NullParameters();
+        warn_initialize_determined = true,
         check_length = true, kwargs...) where {iip}
     if !iscomplete(sys)
         error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEProblem`")
     end
     f, du0, u0, p = process_DEProblem(DAEFunction{iip}, sys, u0map, parammap;
         implicit_dae = true, du0map = du0map, check_length,
-        kwargs...)
+        warn_initialize_determined, kwargs...)
     diffvars = collect_differential_variables(sys)
     sts = unknowns(sys)
     differential_vars = map(Base.Fix2(in, diffvars), sts)
@@ -1496,18 +1511,33 @@ end
 
 function InitializationProblem{iip, specialize}(sys::AbstractODESystem, u0map = [],
         parammap = DiffEqBase.NullParameters();
+        guesses = [],
         check_length = true,
+        warn_initialize_determined = true,
         kwargs...) where {iip, specialize}
     if !iscomplete(sys)
         error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`")
     end
 
-    isys = get_initializesystem(sys)
+    if isempty(u0map)
+        isys = get_initializesystem(sys)
+    else
+        isys = structural_simplify(generate_initializesystem(sys; u0map); fully_determined = false)
+    end
+
     neqs = length(equations(isys))
     nunknown = length(unknowns(isys))
+
+    if warn_initialize_determined && neqs > nunknown
+        @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false."
+    end
+    if warn_initialize_determined && neqs < nunknown
+        @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false."
+    end
+
     if neqs == nunknown
-        NonlinearProblem(isys, u0map, parammap)
+        NonlinearProblem(isys, guesses, parammap)
     else
-        NonlinearLeastSquaresProblem(isys, u0map, parammap)
+        NonlinearLeastSquaresProblem(isys, guesses, parammap)
     end
 end

src/systems/nonlinear/initializesystem.jl

@@ -3,7 +3,9 @@ $(TYPEDSIGNATURES)
 
 Generate `NonlinearSystem` which initializes an ODE problem from specified initial conditions of an `ODESystem`.
 """
-function generate_initializesystem(sys::ODESystem; name = nameof(sys),
+function generate_initializesystem(sys::ODESystem; 
+        u0map = Dict(),
+        name = nameof(sys),
         guesses = Dict(), check_defguess = false, kwargs...)
     sts, eqs = unknowns(sys), equations(sys)
     idxs_diff = isdiffeq.(eqs)
@@ -13,7 +15,7 @@ function generate_initializesystem(sys::ODESystem; name = nameof(sys),
     # Start the equations list with algebraic equations
     eqs_ics = eqs[idxs_alge]
     u0 = Vector{Pair}(undef, 0)
-    defs = defaults(sys)
+    defs = merge(defaults(sys),todict(u0map))
 
     full_states = [sts; getfield.((observed(sys)), :lhs)]
 

src/systems/systemstructure.jl

@@ -615,7 +615,7 @@ function structural_simplify!(state::TearingState, io = nothing; simplify = fals
 end
 
 function _structural_simplify!(state::TearingState, io; simplify = false,
-        check_consistency = true, fully_determined = true, warn_initialize_determined = true,
+        check_consistency = true, fully_determined = true, warn_initialize_determined = false,
         kwargs...)
     check_consistency &= fully_determined
     has_io = io !== nothing
@@ -644,7 +644,7 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
         neqs = length(equations(isys))
         nunknown = length(unknowns(isys))
         if warn_initialize_determined && neqs > nunknown
-            @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares"
+            @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares if $(nunknown - neqs) defaults are not supplied at construction time."
         end
         if warn_initialize_determined && neqs < nunknown
             @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares"

test/initializationsystem.jl

@@ -1,6 +1,49 @@
 using ModelingToolkit, OrdinaryDiffEq, NonlinearSolve, Test
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
+@parameters g
+@variables x(t) y(t) [state_priority = 10] λ(t)
+eqs = [
+       D(D(x)) ~ λ * x
+       D(D(y)) ~ λ * y - g
+       x^2 + y^2 ~ 1
+      ]
+@mtkbuild pend = ODESystem(eqs,t)
+
+initprob = ModelingToolkit.InitializationProblem(pend, [], [g => 1]; guesses = [ModelingToolkit.missing_variable_defaults(pend); x => 1; y => 0.2])
+conditions = getfield.(equations(initprob.f.sys),:rhs)
+
+@test initprob isa NonlinearLeastSquaresProblem
+sol = solve(initprob)
+@test SciMLBase.successful_retcode(sol)
+@test maximum(abs.(sol[conditions])) < 1e-14
+
+initprob = ModelingToolkit.InitializationProblem(pend, [x => 1, y => 0], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend))
+@test initprob isa NonlinearProblem
+sol = solve(initprob)
+@test SciMLBase.successful_retcode(sol)
+@test sol.u == [1.0,0.0,0.0,0.0]
+@test maximum(abs.(sol[conditions])) < 1e-14
+
+initprob = ModelingToolkit.InitializationProblem(pend, [], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend))
+@test initprob isa NonlinearLeastSquaresProblem
+sol = solve(initprob)
+@test !SciMLBase.successful_retcode(sol)
+
+prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], guesses = ModelingToolkit.missing_variable_defaults(pend))
+prob.f.initializeprob isa NonlinearProblem
+sol = solve(prob.f.initializeprob)
+@test maximum(abs.(sol[conditions])) < 1e-14
+sol = solve(prob, Rodas5P())
+@test maximum(abs.(sol[conditions][1])) < 1e-14
+
+prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], guesses = ModelingToolkit.missing_variable_defaults(pend))
+prob.f.initializeprob isa NonlinearLeastSquaresProblem
+sol = solve(prob.f.initializeprob)
+@test maximum(abs.(sol[conditions])) < 1e-14
+sol = solve(prob, Rodas5P())
+@test maximum(abs.(sol[conditions][1])) < 1e-14
+
 @connector Port begin
     p(t)
     dm(t) = 0, [connect = Flow]
@@ -171,34 +214,26 @@ end
 
 @mtkbuild sys = System()
 initprob = ModelingToolkit.InitializationProblem(sys)
+conditions = getfield.(equations(initprob.f.sys),:rhs)
+
 @test initprob isa NonlinearLeastSquaresProblem
 @test length(initprob.u0) == 2
 initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12)
 @test SciMLBase.successful_retcode(initsol)
+@test maximum(abs.(initsol[conditions])) < 1e-14
 
 allinit = unknowns(sys) .=> initsol[unknowns(sys)]
 prob = ODEProblem(sys, allinit, (0, 0.1))
-sol = solve(prob, Rodas5P())
+sol = solve(prob, Rodas5P(), initializealg = BrownFullBasicInit())
 # If initialized incorrectly, then it would be InitialFailure
 @test sol.retcode == SciMLBase.ReturnCode.Unstable
-SciMLBase.has_initializeprob(prob.f)
-
-isys = ModelingToolkit.get_initializesystem(sys)
-unknowns(isys)
-
-initprob = ModelingToolkit.InitializationProblem(sys)
-sol = solve(initprob)
-
-unknowns(sys)
-
-[sys.act.vol₁.dr]
-
-getter = ModelingToolkit.getu(initprob, unknowns(sys)[end-1:end])
-getter(sol)
+@test maximum(abs.(initsol[conditions][1])) < 1e-14
 
-prob.f.initializeprobmap(initsol)
-
-initsol[unknowns(isys)]
+prob = ODEProblem(sys, [], (0, 0.1), check=false)
+sol = solve(prob, Rodas5P())
+# If initialized incorrectly, then it would be InitialFailure
+@test sol.retcode == SciMLBase.ReturnCode.Unstable
+@test maximum(abs.(initsol[conditions][1])) < 1e-14
 
 @connector Flange begin
     dx(t), [guess = 0]
@@ -269,3 +304,28 @@ sol = solve(prob, Rodas5P())
 # If initialized incorrectly, then it would be InitialFailure
 @test sol.retcode == SciMLBase.ReturnCode.Success
 
+prob = ODEProblem(sys, [], (0, 0.1))
+sol = solve(prob, Rodas5P())
+@test sol.retcode == SciMLBase.ReturnCode.Success
+
+### Ensure that non-DAEs still throw for missing variables without the initialize system
+
+@parameters σ ρ β
+@variables x(t) y(t) z(t)
+
+eqs = [D(D(x)) ~ σ * (y - x),
+    D(y) ~ x * (ρ - z) - y,
+    D(z) ~ x * y - β * z]
+
+@mtkbuild sys = ODESystem(eqs, t)
+
+u0 = [D(x) => 2.0,
+    y => 0.0,
+    z => 0.0]
+
+p = [σ => 28.0,
+    ρ => 10.0,
+    β => 8 / 3]
+
+tspan = (0.0, 100.0)
+@test_throws ArgumentError prob = ODEProblem(sys, u0, tspan, p, jac = true)