diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 51a495e4ff..6f0444f2f2 100644 --- a/src/systems/abstractsystem.jl +++ b/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) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index c94166103b..ebea3a82c3 100644 --- a/src/systems/callbacks.jl +++ b/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) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 5f37e524e1..ceaa103288 100644 --- a/test/symbolic_events.jl +++ b/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