simplification of affect systems with observed equations can yield incorrect updates to state variables

[vyudu]

Describe the bug 🐞

Simplifying observed equations in affect systems can lead to unintended explicit updates to state variables when a callback is triggered. All the algebraic equations of the parent system and the observed equations of the parent system are added to the affect system. When the affect system is simplified, some of these observed equations get solved for state variables. Thus the simplified affect system may be fully explicit (i.e. have no actual equations, only observed equations), but this means there will be an explicit update to the state variable that was solved for. Instead of the observed value being updated in response to changes in the state variables, the state variable gets updated according to the previous value of the observed variable (since there are no algebraic equations to be respected).

Expected behavior

Observed equations of the parent system should not be solved for state variables in the affect system.

Minimal Reproducible Example 👇

using ModelingToolkit
t = ModelingToolkit.t_nounits; D = ModelingToolkit.D_nounits
@variables V(t) = -70.0 G(t) = 0.0 jcn(t) [input = true] z(t)
@parameters Eₘ=-70. θ=-50.

ev = (V ≥ θ) => [V ~ Eₘ]

@named lif1 = System([D(V) ~ jcn], t, [V, jcn], [Eₘ, θ]; discrete_events = [ev])
@named lif2 = System([D(V) ~ jcn], t, [V, jcn], [Eₘ, θ]; discrete_events = [ev])
@named qif = System([D(G) ~ G + z, D(z) ~ z, D(V) ~ jcn], t, [V, G, z, jcn], [Eₘ, θ]; discrete_events = [ev])

sys2 = compose(System([lif1.jcn ~ lif2.V + qif.G, lif2.jcn ~ lif1.V + qif.G, qif.jcn ~ lif1.V + lif2.V], t; name = :outer), [lif1, lif2, qif])
sys2 = mtkcompile(sys2)

Then inspecting the update equations for the first discrete event using ModelingToolkit.observed(discrete_events(sys2)[1].affect.system), we get:

4-element Vector{Equation}:
 lif1₊V(t) ~ lif1₊Eₘ
 lif1₊jcn(t) ~ lif2₊jcn(t) - 2lif1₊V(t) + qif₊jcn(t)
 qif₊G(t) ~ lif2₊jcn(t) - lif1₊V(t)
 lif2₊V(t) ~ lif1₊jcn(t) - qif₊G(t)

The observed equation lif2.jcn ~ lif1.V + qif.Ghas been solved for qif.G. Thus the explicit affect will update qif.G using the previous value of lif2.jcn, but what the event should do is change lif1.V and then the value of lif2.jcn should update in response.

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

  [7d9f7c33] Accessors v0.1.42
⌃ [459566f4] DiffEqCallbacks v4.9.0
⌃ [31c24e10] Distributions v0.25.120
  [bcd5d0fe] GraphDynamics v0.4.9 `../GraphDynamics.jl`
  [86223c79] Graphs v1.13.1
  [626554b9] MetaGraphs v0.8.1
  [961ee093] ModelingToolkit v10.26.0
  [c1a47a78] NeurobloxBase v0.1.0 `NeurobloxBase`
  [92933f4c] ProgressMeter v1.11.0
⌃ [731186ca] RecursiveArrayTools v3.37.1
  [189a3867] Reexport v1.2.2
⌃ [0bca4576] SciMLBase v2.120.0
  [2913bbd2] StatsBase v0.34.6
⌅ [d1185830] SymbolicUtils v3.32.0
⌃ [0c5d862f] Symbolics v6.55.0
  [9a3f8284] Random v1.11.0
  [2f01184e] SparseArrays v1.12.0

• Output of versioninfo()

Julia Version 1.12.0
Commit b907bd0600f (2025-10-07 15:42 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 8 × Apple M2
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, apple-m2)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 4 virtual cores)
Environment:
  JULIA_EDITOR = vim
  JULIA_GR_PROVIDER = GR
  JULIA_PKG_PRESERVE_TIERED_INSTALLED = true

Additional context

Add any other context about the problem here.

[AayushSabharwal]

I don't quite follow the chain of reasoning here.

4-element Vector{Equation}:
 lif1₊V(t) ~ lif1₊Eₘ
 lif1₊jcn(t) ~ lif2₊jcn(t) - 2lif1₊V(t) + qif₊jcn(t)
 qif₊G(t) ~ lif2₊jcn(t) - lif1₊V(t)
 lif2₊V(t) ~ lif1₊jcn(t) - qif₊G(t)

These are the observed equations in the affect. There isn't a Pre here, so all these variables should represent the values after the event triggers. As I understand it, you want qif.G to retain its value in the event, and lif2.jcn to be updated as qif.G + lif1.V.

This does not require removing observed equations from the affect. Rather, your affect is underdetermined, and there are multiple solutions that satisfy the constraints you've given. MTK just happens to choose one that you don't want. This is an issue I've been meaning to try and surface, but am not sure how. The problem boils down to the fact that MTK has no idea which variables you want to hold constant, and which to solve for. The solution here is to add qif.G ~ Pre(qif.G) to the affect, which will force it to be held constant and thus solve for lif2.jcn.

[vyudu]

Sorry, I was being unclear, what I meant was there are no algebraic equations in the affect system in this case, so it follows the explicit_affect code generation path, which uses all the previous values on the RHS. So the issue is that turning lif2.jcn ~ lif1.V + qif.G into qif.G ~ lif2.jcn - lif1.V during mtkcompile will generate a different solution - the first updates jcn as we want and the latter does not. The point is taken that this indeterminacy is because the affect is underdetermined, but I think in most cases, the intuitive result is that the observed variable changes in response to changes in the state variable, and is not held constant before and after the affect.

Not sure if adding qif.G ~ Pre(qif.G) to all the systems will work because in Neuroblox it is not determined how many other neurons it will be connected to, and what type they will be, so you can't pass these equations directly to the callback constructor for an individual neuron.qif doesn't exist at the time the event is created for the first two neurons for example.

[AayushSabharwal]

I see. I can agree that in a fully explicit DAE, the obvious behavior is for observed variables to remain observed variables. However, the absence of algebraic equations in the simplified system is not a sufficient condition for excluding observed equations. This also needs to make sure that the user-provided affect equation does not involve non-Pre observed variables. For example, if the affect was V ~ E_m + lif1.jcn then we do need to include the observed equations. This would not be the case for V ~ E_m + Pre(lif1.jcn).