MTK does not allow building initialization systems for un-simplified systems

[SebastianM-C]

Describe the bug 🐞

Let us consider the following system:

using ModelingToolkit

ps = @parameters k = 2 w = 3
sts = @variables x(t) = 1.5

eqs = [
    D(x) ~ -k * x
]
model = structural_simplify(ODESystem(eqs, t, sts, ps; name=:model))

We want to set the value of k to w, change w and then see that k got updated. We don't want changes to the original model, so any model changes have to be out of place.

We can try this in 2 separate ways:

• extending the model with an additional one that adds the new defaults and transforms k in a solvable parameter.
• copy the model and modify the defaults and guesses in-place.

The first option would be:

k′ = Symbolics.wrap(setmetadata(Symbolics.unwrap(k),
    Symbolics.VariableDefaultValue, missing))

extra_sys = ODESystem(Equation[], ModelingToolkit.get_iv(model), [], [k′];
    name = nameof(model),
    description = ModelingToolkit.description(model),
    gui_metadata = ModelingToolkit.get_gui_metadata(model),
    defaults = Dict(k′ => missing),
    guesses = Dict(k′ => 2)
)
Symbolics.metadata(only(parameters(extra_sys)))
model′ = complete(extend(extra_sys, model))

defaults(model′)[k′] # missing
Symbolics.getmetadata(Symbolics.value(k′), Symbolics.VariableDefaultValue) # missing
Symbolics.getmetadata(Symbolics.value(model′.k), Symbolics.VariableDefaultValue) # 2

prob = ODEProblem(model′, [], (0.,1), [k=>w])
prob.f.initialization_data # nothing
(u0, p) = setsym_oop(prob, w)(prob, 3.1)
prob_new = remake(prob; u0, p)
init(prob_new).ps[k] # 3.0

Note that we have a couple of strange things happening:

• the new default that we introduce is with a variable where we explicitly set the metadata, but when we look at the k inside the system, we see that we have the old value
• no initialization problem is built, even if we have a parameter that should satisfy the conditions

This results in the parameter update being ignored.

Now, for the second, we do it via via deepcopy:

model2 = deepcopy(model)

guesses(model2)[k] = 2
defaults(model2)[k] = missing

prob2 = ODEProblem(model2, [], (0., 1), [k => w])
prob2.f.initialization_data # SciMLBase.OverrideInitData{NonlinearProblem{...}}
(u0, p) = setsym_oop(prob2, w)(prob2, 3.1)
prob2_new = remake(prob2; u0, p)
init(prob2_new).ps[k] # 3.1

While this works as expected, the issue is that deepcopy sometimes triggers initialization issues for some reason, leading to #3307

Expected behavior

A clear and concise description of what you expected to happen.

Minimal Reproducible Example 👇

using ModelingToolkit, Symbolics
using OrdinaryDiffEqDefault

ps = @parameters k = 2 w = 3
sts = @variables x(t) = 1.5

eqs = [
    D(x) ~ -k * x
]
model = structural_simplify(ODESystem(eqs, t, sts, ps; name=:model))

k′ = Symbolics.wrap(setmetadata(Symbolics.unwrap(k),
    Symbolics.VariableDefaultValue, missing))

extra_sys = ODESystem(Equation[], ModelingToolkit.get_iv(model), [], [k′];
    name = nameof(model),
    description = ModelingToolkit.description(model),
    gui_metadata = ModelingToolkit.get_gui_metadata(model),
    defaults = Dict(k′ => missing),
    guesses = Dict(k′ => 2)
)
Symbolics.metadata(only(parameters(extra_sys)))
model′ = complete(extend(extra_sys, model))

defaults(model′)[k′] # missing
Symbolics.getmetadata(Symbolics.value(k′), Symbolics.VariableDefaultValue) # missing
Symbolics.getmetadata(Symbolics.value(model′.k), Symbolics.VariableDefaultValue) # 2

prob = ODEProblem(model′, [], (0.,1), [k=>w])
prob.f.initialization_data # nothing
(u0, p) = setsym_oop(prob, w)(prob, 3.1)
prob_new = remake(prob; u0, p)
init(prob_new).ps[k] # 3.0

Error & Stacktrace ⚠️

There is no error, but the results are wrong since the parameter k is not solved for during initialization.

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

(dev) pkg> st -m ModelingToolkit Symbolics SymbolicUtils
Status `~/dev/Manifest-v1.11.toml`
  [961ee093] ModelingToolkit v9.60.0
  [d1185830] SymbolicUtils v3.8.1
  [0c5d862f] Symbolics v6.22.1

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

• Output of versioninfo()

julia> versioninfo()
Julia Version 1.11.2
Commit 5e9a32e7af2 (2024-12-01 20:02 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × Intel(R) Core(TM) i9-14900K
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, alderlake)
Threads: 32 default, 0 interactive, 16 GC (on 32 virtual cores)
Environment:
  JULIA_EDITOR = code

Additional context

A workaround to this issue is to use deepcopy, but that hits #3307.

[AayushSabharwal]

That's not how this should be done. Defaults in variable metadata are just a convenience feature. The source of truth for defaults in a completed system is defaults(sys). If you wish to update the defaults in an out of place manner,

newsys = ODESystem(Equation[], t, [], [k]; defaults = [k => missing], guesses = [k => 2], name = :foo)
model2 = extend(newsys, model) # `newsys` needs to be first

[SebastianM-C]

I initially tried that, but I saw that no initialization problem was created and this is what prompted me to look at metadata and try to force a change there, but that didn't work either.

Note that in my example above I updated both the defaults (and guesses) and the metadata. In both cases (the code below and my above example), the defaults for the updated model looks correct, but the metadata of the changed parameter does not coincide with what defaults says. I'm not sure if that is in any way related to the initialization problem not being created.

using ModelingToolkit
using ModelingToolkit: D_nounits as D, t_nounits as t
using SymbolicIndexingInterface
using SciMLBase, OrdinaryDiffEqDefault

ps = @parameters k = 2 w = 3
sts = @variables x(t) = 1.5

eqs = [
    D(x) ~ -k * x
]
model = structural_simplify(ODESystem(eqs, t, sts, ps; name=:model))

newsys = ODESystem(Equation[], t, [], [k]; defaults=[k => missing], guesses=[k => 2], name=:foo)
model2 = complete(extend(newsys, model)) # `newsys` needs to be first

prob = ODEProblem(model2, [], (0., 1), [k => w])
prob.f.initialization_data # nothing
(u0, p) = setsym_oop(prob, w)(prob, 3.1)
prob_new = remake(prob; u0, p)
init(prob_new).ps[k] # 3.0

[AayushSabharwal]

Initialization not being created is a different issue. That happens because according to MTK model2 is not structural_simplified since MTK.get_tearing_state(model2) === nothing. I'm working on addressing it, but the issue is more annoying than it first appears to be.

[AayushSabharwal]

As a workaround, @set! model2.tearing_state = MTK.get_tearing_state(model) should make MTK happy

[SebastianM-C]

Thanks! Should I close this issue then?

[AayushSabharwal]

I think we can keep it up as a tracking issue for allowing initialization without structural_simplify. I don't think we have one yet.

[SebastianM-C]

I see that on a hierarchical model the preamble code in the generated function (where the observed are defined) gets dropped when I apply this transformation. Am I missing something?