fix: only add observed equations to affect system if there are algebraic equations in the system

src/systems/abstractsystem.jl

@@ -645,15 +645,17 @@ function complete(
         if add_initial_parameters
             sys = add_initialization_parameters(sys; split)
         end
+        alg_eqs = alg_equations(sys)
+        obs_eqs = observed(sys)
         if has_continuous_events(sys) && is_time_dependent(sys)
             @set! sys.continuous_events = complete.(
                 get_continuous_events(sys); iv = get_iv(sys),
-                alg_eqs = [alg_equations(sys); observed(sys)])
+                alg_eqs, obs_eqs)
         end
         if has_discrete_events(sys) && is_time_dependent(sys)
             @set! sys.discrete_events = complete.(
                 get_discrete_events(sys); iv = get_iv(sys),
-                alg_eqs = [alg_equations(sys); observed(sys)])
+                alg_eqs, obs_eqs)
         end
     end
     if split && has_index_cache(sys)

src/systems/callbacks.jl

@@ -61,7 +61,7 @@ function AffectSystem(spec::SymbolicAffect; iv = nothing, alg_eqs = Equation[],
 end
 
 function AffectSystem(affect::Vector{Equation}; discrete_parameters = Any[],
-        iv = nothing, alg_eqs::Vector{Equation} = Equation[], warn_no_algebraic = true, kwargs...)
+        iv = nothing, alg_eqs::Vector{Equation} = Equation[], obs_eqs::Vector{Equation} = Equation[], warn_no_algebraic = true, kwargs...)
     isempty(affect) && return nothing
     if isnothing(iv)
         iv = t_nounits
@@ -70,6 +70,7 @@ function AffectSystem(affect::Vector{Equation}; discrete_parameters = Any[],
 
     discrete_parameters isa AbstractVector || (discrete_parameters = [discrete_parameters])
     discrete_parameters = unwrap.(discrete_parameters)
+    obs_vars = [unwrap(eq.lhs) for eq in obs_eqs]
 
     for p in discrete_parameters
         occursin(unwrap(iv), unwrap(p)) ||
@@ -79,6 +80,7 @@ function AffectSystem(affect::Vector{Equation}; discrete_parameters = Any[],
     dvs = OrderedSet()
     params = OrderedSet()
     _varsbuf = Set()
+    eqs_contain_obs = false
     for eq in affect
         if !haspre(eq) && !(symbolic_type(eq.rhs) === NotSymbolic() ||
              symbolic_type(eq.lhs) === NotSymbolic())
@@ -91,6 +93,16 @@ function AffectSystem(affect::Vector{Equation}; discrete_parameters = Any[],
         union!(params, _varsbuf)
         diffvs = collect_applied_operators(eq, Differential)
         union!(dvs, diffvs)
+
+        for var in obs_vars
+            if occursin(var, unwrap(eq.rhs))
+                eqs_contain_obs = true
+            end
+        end
+    end
+
+    if eqs_contain_obs
+        alg_eqs = [alg_eqs; obs_eqs]
     end
     for eq in alg_eqs
         collect_vars!(dvs, params, eq, iv)

test/symbolic_events.jl

@@ -1440,3 +1440,46 @@ end
     @mtkcompile sys = MWE()
     @test_nowarn ODEProblem(sys, [], (0.0, 1.0))
 end
+
+@testset "Test explicit updates correctly handle observed variables." begin
+    t = ModelingToolkit.t_nounits; D = ModelingToolkit.D_nounits
+    @variables V(t)=-70.0 G(t)=0.0 jcn(t) [input = true] z(t)= 0.0
+    @parameters Eₘ=-70.0 θ=-50.0
+
+    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)
+
+    ev = discrete_events(sys2)[1]
+    @test length(observed(ev.affect.system)) == 1
+    prob = ODEProblem(sys2, [], (0., 2.))
+    sol = solve(prob)
+
+    first_ev_idx = findfirst(i -> sol.t[i] == sol.t[i+1], 1:length(sol.t))
+    t_e = sol.t[first_ev_idx]
+    @test sol(t_e - eps(), idxs = [lif2.V, lif1.V]) ≈ sol(t_e + eps(), idxs = [lif2.V, lif1.V])
+
+    # Add observed equations if explicit observed
+    ev = (V ≥ θ) => [V ~ -jcn]
+    @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])
+
+    sys3 = 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])
+    sys3 = mtkcompile(sys3)
+
+    ev = discrete_events(sys3)[1]
+    @test length(observed(ev.affect.system)) == 4
+    prob = ODEProblem(sys3, [], (0., 2.))
+    sol = solve(prob)
+
+    first_ev_idx = findfirst(i -> sol.t[i] == sol.t[i+1], 1:length(sol.t))
+    t_e = sol.t[first_ev_idx]
+    @test sol(t_e + eps(), idxs = qif.V) ≈ -sol(t_e + eps(), idxs = qif.jcn)
+    @test sol(t_e - eps(), idxs = [lif1.V, lif2.V]) ≈ sol(t_e + eps(), idxs = [lif1.V, lif2.V])
+end