Providing parameter override errors with "Initialization expression currently not supported"

[baggepinnen]

The following (minimized) example from https://help.juliahub.com/juliasimcontrol/dev/robust_control/ has been broken on some recent release of MTK, the error is

ERROR: Initialization expression sys₊k is currently not supported. If its a higher order derivative expression, then only the dummy derivative expressions are supported.

using ModelingToolkit, MonteCarloMeasurements, OrdinaryDiffEq

T = typeof(1 ± 1)
@parameters t k::T d
@variables x(t)::T v(t)::T
D = Differential(t)
eqs = [D(x) ~ v, D(v) ~ -k*x - d*v - 9.82]
@named sys = ODESystem(eqs, t)
prob = ODEProblem(complete(sys), [x => 0.0 ± 0.1, v => 0.0, k => 10 ± 1, d => 1], (0.0, 10.0))
sol = solve(prob, Tsit5())

T = typeof(-1..1)
@parameters Δ::T=(-1..1) # The uncertain element is between -1 and 1
@variables y(t)=0
connections = [
    y ~ sys.x + Δ
]

@named uncertain_sys = ODESystem(connections, t, systems=[sys])
initial_condition = [sys.x => zero(T), sys.v => 0.0, sys.k => 10, sys.d => 1]
prob = ODEProblem(structural_simplify(uncertain_sys, split=false), initial_condition, (0.0, 10.0)) # This line errors
sol = solve(prob, Tsit5())

here's a failing build log
https://github.com/JuliaComputing/JuliaSimControl.jl/actions/runs/11119591661/job/30899449776#step:6:442

[AayushSabharwal]

This is because we didn't build initialization systems for defaults of observed variables before, but now we do. Since y has a default, the initialization system is built. It errors because u0map contains parameters.

[ChrisRackauckas]

Yes, it was simply wrong before #3032. That PR cleaned up a lot of cases that should use the initialization that did not before. Now that's the case, it errors because you have an issue with the u0map is not correctly built.

[AayushSabharwal]

There are several layers here. Fixing the u0map reveals that our calculation of resid_prototype for the NonlinearFunction is problematic. Since u0 === nothing (the NonlinearSystem has 0 unknowns and 1 equation) it takes resid_prototype = zeros(Float64, 1). However, the parameters are Particles{Float64, 2000} from MonteCarloMeasurements. Thus, it errors when trying to solve the NonlinearProblem since the in-place function tries to set a Vector{Float64} with a Particles{Float64, 2000}. Fixing this leads to an InitialFailure, because given these observed equations:

 sys₊x(t) ~ 0.0
 sys₊v(t) ~ 0.0
 y(t) ~ 0.0

and this equation:

0 ~ sys₊x(t) - y(t) + Δ

and Δ = -1.6098200000000002e-18 ± 0.577, norm(resid) is 0.5 ± 0.289 which is above the abstol of 1e-6

So apart from fixing NonlinearFunction, y shouldn't have a default.

[ChrisRackauckas]

Yes, there's a user error there for sure. Though we should do a type promotion on the resid_prototype similar to the dual handling based on the u0 and tunables type.

[ChrisRackauckas]

@AayushSabharwal then anything else here?

[AayushSabharwal]

No