How to prevent breakage from changes to `structural_simplify`

[baggepinnen]

I very often run into problems where tests or documentation builds suddenly start failing because a change has been made somewhere in the simplification pipeline that changes what variables I need to provide initial conditions for. As an example, I just started getting the error

SymbolicUtils.Symbolic{Real}[link1₊x1ˍtt(t), link1₊y_cmˍtt(t)] are missing from the variable map.

in a documentation example that has previously been running just fine. I'm now being requested to provide initial conditions for the second derivatives of some of the variables in my system. Technically, any change to structural_simplify that would force the user to provide such initial conditions where they were not previously required is a breaking change. Is there a good way for the user to guard themselves from such changes to not have tests broken every now and then?

[ChrisRackauckas]

Make your tests be downstream tests to MTK?

Technically, any change to structural_simplify that would force the user to provide such initial conditions where they were not previously required is a breaking change.

In theory, you only need the minimal set, and it should be able to deduce some of them from the others. In practice that part isn't implemented yet so it's just waiting in the wings.

[finmod]

Along the same issue, this DAE model of constrained dynamics used to work before Julia v.1.8.4 and some modfifications in MTK in December. Now, after structural_simplify it breaks down on states[A1_ode] and later using ODAEProblem it gives the error mentioned by @baggepinnen

SymbolicUtils.Symbolic{Real}[link1₊x1ˍtt(t), link1₊y_cmˍtt(t)] are missing from the variable map.

using DifferentialEquations
using OrdinaryDiffEq
using ModelingToolkit
using Plots; gr()

ps = @parameters α β dp̂ γ K₀ l₀ L̂ Mᶠ₀ Mʰ₀ M̂ʰ p₀ S₀ Ŝ w₀ μ₁ μ₂ μ₃ μ₄ μ₅ μ₆ μ₇ μ₈ μ₉ μ₁₀ 𝛌₁ 𝛌₂ 𝛌₃

@variables t 
states = @variables  c(t) K(t) l(t) Mᶠ(t) Mʰ(t) p(t) S(t) w(t) k(t) mᶠ(t) mʰ(t) s(t)
@variables dep(t) y(t) Uᶠ(t) Uʰ(t) Zₕ(t) Zₚ(t) Zᵣ(t) 𝛌₁(t) 𝛌₂(t) 𝛌₃(t)

Dₜ = Differential(t) 
                  
dep = dp̂*K
y   = β*K^(1 - α) * l^(α)

Uᶠ  = p*y - w*l - (Ŝ - S)^2 
Uʰ  = c^(γ) - (L̂ - l)^2 - (M̂ʰ - Mʰ)^2 

Zₕ  = w*l - p*c - mʰ
Zₚ  = p*c - w*l - mᶠ            
Zᵣ  = y - c - k - s - dep


DLP = [Uᶠ, Uʰ, 𝛌₁*Zₕ, 𝛌₂*Zₚ, 𝛌₃*Zᵣ]

∂DLPᵀ = transpose(Symbolics.jacobian(DLP, states))


eqs = [Dₜ(c)  ~                 μ₁*∂DLPᵀ[1,2] + ∂DLPᵀ[1,3] + ∂DLPᵀ[1,4] + ∂DLPᵀ[1,5],
       Dₜ(K)  ~ μ₂*∂DLPᵀ[2,1]                 + ∂DLPᵀ[2,3] + ∂DLPᵀ[2,4] + ∂DLPᵀ[9,5],
       Dₜ(l)  ~ μ₃*∂DLPᵀ[3,1] + μ₄*∂DLPᵀ[3,2] + ∂DLPᵀ[3,3] + ∂DLPᵀ[3,4] + ∂DLPᵀ[3,5],
       Dₜ(Mᶠ) ~ μ₅*∂DLPᵀ[4,1]                 + ∂DLPᵀ[4,3] + ∂DLPᵀ[10,4]+ ∂DLPᵀ[4,5],
       Dₜ(Mʰ) ~                 μ₆*∂DLPᵀ[5,2] + ∂DLPᵀ[11,3]+ ∂DLPᵀ[5,4] + ∂DLPᵀ[5,5],
       Dₜ(p)  ~ μ₇*∂DLPᵀ[6,1]                 + ∂DLPᵀ[6,3] + ∂DLPᵀ[6,4] + ∂DLPᵀ[6,5],
       Dₜ(S)  ~ μ₈*∂DLPᵀ[7,1]                 + ∂DLPᵀ[7,3] + ∂DLPᵀ[7,4] + ∂DLPᵀ[12,5],
       Dₜ(w)  ~ μ₉*∂DLPᵀ[8,1]                 + ∂DLPᵀ[8,3] + ∂DLPᵀ[8,4] + ∂DLPᵀ[8,5],
       Dₜ(K)  ~ k,                                
       Dₜ(Mᶠ) ~ mᶠ,
       Dₜ(Mʰ) ~ mʰ,
       Dₜ(S)  ~ s,
       0.     ~ Zₕ,
       0.     ~ Zₚ,
       0.     ~ Zᵣ]

defs = Dict(Uᶠ => 1.0, Uʰ => 1.0, dep => 1.0, y => 1.0, c => β*K₀^(1 - α) * l₀^(α), K => K₀, l => l₀, Mᶠ => Mᶠ₀, Mʰ => Mʰ₀, p => p₀, S => S₀, w => w₀, k => 0.1, mᶠ => 0, mʰ => 0, s => 0, 𝛌₁ => 0.01, 𝛌₂ => 0.01, 𝛌₃ => 0.02,
       α => 0.5, β => 1.0, dp̂ => 0.1, γ => 0.5, K₀ => 1.6, l₀ => 1.0, L̂ => 1.0, Mᶠ₀ => 1.05, Mʰ₀ => 1.2, M̂ʰ => 1.0, p₀ => 0.8, S₀ => 0.4, Ŝ  => 0.1, w₀ => 1.2, μ₁ => 1.0, μ₂ => 1.0, μ₃ => 1.0, μ₄ => 1.0, μ₅ => 1.0, μ₆ => 1.0, μ₇ => 1.0, μ₈ => 1.0, μ₉ => 0.0, μ₁₀ => 0.0)

@named A1_dae = ModelingToolkit.ODESystem(eqs, t, states, ps; defaults = defs, checks = false)

A1_ode = structural_simplify(dae_index_lowering(A1_dae))
parameters(A1_ode)
observed(A1_ode)
states[A1_ode]
```Julia

[ChrisRackauckas]

Yes, the real issue is that we're missing an initialization module which handles the issue that defining an observable's initial condition gives constraints on the initial condition of the solver states. Pointing to structural simplify changes is an X-Y problem: this just needs the initialization to exist.

[baggepinnen]

Just providing better initialization is not quite sufficient, some applications call for a user-chosen state realization for reasons other than initialization. Optimal control is one such application, where the user might want a particular realization for numerical reasons etc. #2048 is an attempt at providing such options to the user.