Callbacks break systems with initialization

[baggepinnen]

Firing of a callback resets the state to the initial state

[baggepinnen]

MWE

Three ingredients are needed to trigger the reset:

1. a mass-matrix DAE (an algebraic constraint), so the integrator reinitializes after a continuous callback;
2. a non-trivial / implicit initialization (initialization_eqs + guesses), so the problem carries a real initialization system (prob.f.initialization_data !== nothing);
3. a plain SciMLBase.ContinuousCallback with a no-op affect (a user-supplied callback, not MTK's SymbolicContinuousCallback).

using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D
const SII = ModelingToolkit.SymbolicIndexingInterface

@variables x(t) y(t)
eqs = [D(x) ~ y
       0 ~ y^2 - x]                          # algebraic constraint -> DAE
@mtkcompile sys = System(eqs, t; initialization_eqs = [x ~ 1.0])

prob = ODEProblem(sys, [], (0.0, 8.0); guesses = [x => 1.0, y => 1.0])

# A plain ContinuousCallback whose affect does NOTHING -- it only places a solver
# node when x crosses 3. It changes no state and does not call `u_modified!`.
get_x = SII.getu(sys, x)
cond(u, t, integ) = get_x(SII.ProblemState(; u = u, p = integ.p, t = t)) - 3.0
cb = ContinuousCallback(cond, integ -> nothing)

sol    = solve(prob, Rodas5P())
sol_cb = solve(prob, Rodas5P(); callback = cb)

@show sol(4.0;    idxs = x)   # 9.0  (correct)
@show sol_cb(4.0; idxs = x)   # ~2.36

Observed vs expected

x should integrate freely as D(x) = y = sqrt(x), reaching x(4) ≈ 9. The no-op callback should not change the trajectory. Instead, every time the callback fires (x crosses 3), x is reset toward its initial value x = 1 and never gets past 3:

x(0)            = 1.0
x(4) no cb      = 9.0     # correct
x(4) with no-op = 2.36    # reset toward the initial state

Root cause

The problem carries an OverrideInitData, and MTK sets the integrator's initializealg to OverrideInit. When the plain ContinuousCallback fires, the DAE integrator reinitializes and re-runs the full initialization problem, which re-applies initialization_eqs (here x ~ 1) to the differential variable, snapping it back to the initial value.

MTK's own SymbolicContinuousCallback does not hit this because it sets reinitializealg = NoInit(). A user-supplied raw ContinuousCallback has no such protection, so effectively any callback resets the state. (My real model uses such a callback only to put a solver node at an ifelse kink — no state change, no terminate! — and the state still resets.)

Versions

Reproduced on MTK 11.26.8 / ModelingToolkitBase 1.42.2 and current master MTK 11.28.2 / ModelingToolkitBase 1.45.0, with both Rodas5P() and FBDF().

[ChrisRackauckas]

Yes it needs a different initialization for the callbacks that drops the things that are only about the initial point.