Discrete event on a system with algebraic equations changes a differential state instead of the algebraic variable #5200

Summary

A SymbolicDiscreteCallback whose affect changes only a discrete parameter also changes a differential state when the system has an algebraic relation between that state and an algebraic variable. AffectSystem appends the parent's algebraic and observed equations to the affect and then runs mtkcompile on the underdetermined affect system (ModelingToolkitBase/src/systems/callbacks.jl, AffectSystem(affect::Vector{Equation}; ...)). That compile solves 0 ~ 2λ + v - 2on for v, holding λ at its pre-event value, when it should hold v and recompute λ. The integration then continues from a wrong state, with no warning.

Reproducer

using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D
@variables x(t) v(t) λ(t) [guess = 0.0]
@discretes on(t) = 1.0
eqs = [D(x) ~ v, D(v) ~ -λ, 0 ~ 2λ + v - 2on]
ev = ModelingToolkit.SymbolicDiscreteCallback([0.5], [on ~ 0.0]; discrete_parameters = [on])
@mtkcompile sys = System(eqs, t, [x, v, λ], [on]; discrete_events = [ev])
@show unknowns(sys)
prob = ODEProblem(sys, [x => 0.0, v => 1.0], (0.0, 1.0))
sol = solve(prob, Rodas5P(); abstol = 1e-10, reltol = 1e-10)
i = findall(==(0.5), sol.t)
println("state at t = 0.5 before / after the event: ", sol.u[i[1]], " / ", sol.u[i[2]])

Output:

unknowns(sys) = [x(t), v(t)]
state at t = 0.5 before / after the event: [0.4319491666245227, 0.7159745833122585] / [0.4319491666245227, -1.2840254166877414]

v is a differential state, so it should stay at 0.71597 across the event, and λ = (2on - v)/2 should jump from 0.642 to -0.358. What happens is that λ stays at 0.642 and v is set to 2·0 - 2·0.642 = -1.284.

The same happens when λ stays an unknown (a mass-matrix DAE), for example in the damped mass-spring system of SciMLBenchmarks DAE/LinearDAE.jmd, with unknowns [x1, v1, λ]. After an event that switches the input, v1 changes from 1.00245 to -0.98542 and λ keeps its old value.

Second problem: out-of-place SVector problems

With the same system built as ODEProblem{false}(sys, SA[x => 0.0, v => 1.0], (0.0, 1.0)), the event throws:

setindex!(::SVector{2, Float64}, value, ::Int) is not defined.

from ExplicitAffect → u_up!(integ). The generated update writes integ.u in place.

Related

#4661 covers a different reinitialization problem. A plain DiscreteCallback that sets integ.ps[on] = 0.0 on this kind of system, with the default reinitialization, resets x and v to their t = 0 values. I reproduced that on the versions below too. With initializealg = BrownFullBasicInit() on the callback, the in-place problem behaves correctly.

Versions

ModelingToolkit v11.45.1, ModelingToolkitBase v1.77.0, OrdinaryDiffEq v7.8.1 (the latest registered as of 2026-09-23). The same results come from ModelingToolkit v11.43.1 / ModelingToolkitBase v1.71.2.

Julia Version 1.11.9
Commit 53a02c0720c (2026-02-06 00:27 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 128 × AMD EPYC 7502 32-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, znver2)
Threads: 1 default, 0 interactive, 1 GC (on 128 virtual cores)

🤖 Filed by Claude Code (model: claude-opus-5-5[1m]) — https://claude.ai/code/session_01Vrsu4PESdABpNaxYiTBVfd