Merge pull request #2152 from SciML/bgc/dummy_init

missing variable defaults function

Author: YingboMa

docs/src/basics/FAQ.md

@@ -96,3 +96,36 @@ function loss(p)
     sum(abs2, sol)
 end
 ```
+
+# ERROR: ArgumentError: SymbolicUtils.BasicSymbolic{Real}[xˍt(t)] are missing from the variable map.
+
+This error can come up after running `structural_simplify` on a system that generates dummy derivatives (i.e. variables with `ˍt`).  For example, here even though all the variables are defined with initial values, the `ODEProblem` generation will throw an error that defaults are missing from the variable map.
+
+```
+@variables t
+sts = @variables x1(t)=0.0 x2(t)=0.0 x3(t)=0.0 x4(t)=0.0
+D = Differential(t)
+eqs = [x1 + x2 + 1 ~ 0
+       x1 + x2 + x3 + 2 ~ 0
+       x1 + D(x3) + x4 + 3 ~ 0
+       2 * D(D(x1)) + D(D(x2)) + D(D(x3)) + D(x4) + 4 ~ 0]
+@named sys = ODESystem(eqs, t)
+sys = structural_simplify(sys)
+prob = ODEProblem(sys, [], (0,1))
+```
+
+We can solve this problem by using the `missing_variable_defaults()` function
+
+```
+prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0,1))
+```
+
+This function provides 0 for the default values, which is a safe assumption for dummy derivatives of most models.  However, the 2nd arument allows for a different default value or values to be used if needed.
+
+```
+julia> ModelingToolkit.missing_variable_defaults(sys, [1,2,3])
+3-element Vector{Pair}:
+  x1ˍt(t) => 1
+ x2ˍtt(t) => 2
+ x3ˍtt(t) => 3
+```

src/systems/abstractsystem.jl

@@ -1756,3 +1756,31 @@ Base.:(∘)(sys1::AbstractSystem, sys2::AbstractSystem) = compose(sys1, sys2)
 function UnPack.unpack(sys::ModelingToolkit.AbstractSystem, ::Val{p}) where {p}
     getproperty(sys, p; namespace = false)
 end
+
+"""
+    missing_variable_defaults(sys::AbstractSystem, default = 0.0)
+
+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)
+    varmap = get_defaults(sys)
+    varmap = Dict(Symbolics.diff2term(value(k)) => value(varmap[k]) for k in keys(varmap))
+    missingvars = setdiff(states(sys), keys(varmap))
+    ds = Pair[]
+
+    n = length(missingvars)
+
+    if default isa Vector
+        @assert length(default)==n "`default` size ($(length(default))) should match the number of missing variables: $n"
+    end
+
+    for (i, missingvar) in enumerate(missingvars)
+        if default isa Vector
+            push!(ds, missingvar => default[i])
+        else
+            push!(ds, missingvar => default)
+        end
+    end
+
+    return ds
+end

test/odesystem.jl

@@ -853,6 +853,18 @@ let
 
     @named sys = ODESystem(eqs, t, u, ps)
     @test_nowarn simpsys = structural_simplify(sys)
+
+    sys = structural_simplify(sys)
+
+    u0 = ModelingToolkit.missing_variable_defaults(sys)
+    u0_expected = Pair[s => 0.0 for s in states(sys)]
+    @test string(u0) == string(u0_expected)
+
+    u0 = ModelingToolkit.missing_variable_defaults(sys, [1, 2])
+    u0_expected = Pair[s => i for (i, s) in enumerate(states(sys))]
+    @test string(u0) == string(u0_expected)
+
+    @test_nowarn ODEProblem(sys, u0, (0, 1))
 end
 
 # https://github.com/SciML/ModelingToolkit.jl/issues/1583