Continuous Callback fails on observed variable

[dodoplus]

Take a look at this slightly modified code from the tutorial:

using ModelingToolkit, Plots, DifferentialEquations

@variables t
@connector function Pin(;name)
    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
    ODESystem(Equation[], t, sts, []; name=name)
end

function Ground(;name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name=name), g)
end

function OnePort(;name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

function Resistor(;name, R = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R=R
    eqs = [
           v ~ i * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters C=C
    D = Differential(t)
    eqs = [
           D(v) ~ i / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name, continuous_events = [i ~ 0] => [C ~ C-v]), oneport)
end

function ConstantVoltage(;name, V = 1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V=V
    eqs = [
           V ~ v
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

R = 1.0
C = 1.0
V = 1.0
@named resistor = Resistor(R=R)
@named capacitor = Capacitor(C=C)
@named source = ConstantVoltage(V=V)
@named ground = Ground()

rc_eqs = [
          connect(source.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source.n)
          connect(capacitor.n, ground.g)
         ]

@named _rc_model = ODESystem(rc_eqs, t)
@named rc_model = compose(_rc_model,
                          [resistor, capacitor, source, ground])
sys = structural_simplify(rc_model)
u0 = [
      capacitor.v => 0.0
     ]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())
plot(sol)

The only difference is the added continuous_events = [i ~ 0] => [C ~ C-v] in the capacitor code.

However, i is eliminated in structural_simplify:

states(sys)
# 1-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
# capacitor₊v(t)

Running the code above therefore throws an exception:

UndefVarError: capacitor₊i not defined

Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/SymbolicUtils/9IZDP/src/code.jl:394 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/Symbolics/8NTj0/src/build_function.jl:469 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/SymbolicUtils/9IZDP/src/code.jl:351 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:129 [inlined]
  [5] macro expansion
...

Note that the same code with v instead of i works fine.

[baggepinnen]

#1396

one solution is to prevent variables that appear in callbacks to be moved to observerd

[dodoplus]

Is there a workaround?

I think that designating variables as "protected" from elimination is a good idea.

[YingboMa]

You can always mark a variable as irreducible, e.g. @variables x [irreducible=true]

[ValentinKaisermayer]

The value could also be computed in the callback via prob.f.observed. Could this be done automatically if an observed variable is detected in the callback? I feel that such implementation details should be handled automatically or otherwise similar questions will come up very often.

[wsphillips]

I've been looking at this. I'll leave a couple notes for possible correction/feedback from @YingboMa and/or to potentially help anyone else that attempts to fix this:

The symbolics used by callbacks (both conditions and affects) probably need to have a preface and substitutions passed down to Symbolics.build_function similar to the main codegen paths. However there's two things that complicate blindly reapplying the ones made for the main system of equations:

1. There may be observables/substitutions that are excluded and/or unnecessary for the callback function bodies. This means either some get left out, extra observables are being inefficiently recomputed, or both.

The value could also be computed in the callback via prob.f.observed. Could this be done automatically if an observed variable is detected in the callback?

We might be able to pregenerate the code needed with some clever reuse. We can take advantage of the fact that prob.f.observed comes wrapped in a let block with all the values needed for codegen bundled...

ModelingToolkit.jl/src/structural_transformation/codegen.jl

Lines 313 to 330 in a09b914

	observedfun = let state = state,
	dict = Dict(),
	is_solver_state_idxs = insorted.(1:length(fullvars), (states_idxs,)),
	assignments = assignments,
	deps = (deps, invdeps),
	sol_states = sol_states,
	var2assignment = var2assignment
	function generated_observed(obsvar, u, p, t)
	obs = get!(dict, value(obsvar)) do
	build_observed_function(state, obsvar, var_eq_matching, var_sccs,
	is_solver_state_idxs, assignments, deps,
	sol_states, var2assignment,
	checkbounds = checkbounds)
	end
	obs(u, p, t)
	end
	end

Both build_observed_function and build_explicit_observed_function have optional keywords to return an expression instead of a RGF, so they could be called using the fields embedded in prob.f.observed to generate the expression blocks needed for solving arbitrary observed states.

2. Discrete callback conditions are of type Num. To reuse Symbolics.build_function we likely want to take advantage of the keyword arguments states and postprocess_fbody. Right now, those are only supported for vectors of RHS symbolics. I didn't look at that long enough to determine whether it's an easy extension for Symbolics to do or not.