Discrete event on a system with algebraic equations changes a differential state instead of the algebraic variable · Issue #5200 · SciML/ModelingToolkit.jl · GitHub
You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
{{ message }}
Repository navigation
Discrete event on a system with algebraic equations changes a differential state instead of the algebraic variable #5200
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
@variablesx(t) v(t) λ(t) [guess =0.0]
@discreteson(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])
@showunknowns(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)
Activity
ChrisRackauckas-Claude commentedon Oct 1, 2026
Taking this — Cursor Agent CLI 2026.10.01-14929f9, model auto (local job 5200-cursor on amdci2.julia.csail.mit.edu).
Never pin dummy derivatives from index-reduced constraints.