Crossing conditions in ModelingToolkit @mtkmodel macro

Specific Domains Modelling & Simulations
,
1

I want to add a continuous event trigger to my @mtkmodel, with a downcrossing condition. Without the crossing direction condition, I think this should look like

@mtkmodel odesys begin
    @variables begin
        X(t)
        Y(t)
    end
...
    @continuous_events begin
        [X ~ X_thresh] => [Y ~ Y/2]
    end
end

Is there a way to specify downcrossings within the macro, or do I have to abandon @mtkmodel and construct everything explicitly/procedurally by building an ODESystem passing in arrays of equations, variables, events, etc.?

2

I haven’t found any way of doing this automatically within the macro system; I think the macros automatically guarantee that all continuous events trigger on both up- and down-crossings.

The next-best solution I found was to pull info out of the @mtkmodel system with the incorrect positive+negative affects, manually construct a downcrossing event with only a negative affect, and then manually rebuild the ODE system with it. I also had to hack around with some equation substitutions to handle the external forcing code my model has; the substitute code below may not be necessary for all models.

Code:

@mtkbuild odesys_ref = odesys()

# create a continuous event callback that triggers on the old model's threshold crossing but only takes action on downcrossings
ce = continuous_events(odesys_ref)[1] # old model's continuous event definition
cb_down = ModelingToolkit.SymbolicContinuousCallback(
    eqs = ce.eqs,
    affect = Equation[],           # positive-edge affect = nothing (was ce.affect in old model); turn off upcrossing trigger
    affect_neg = ce.affect,       # negative-edge affect (ce.affect); downcrossing trigger is what the old model used
    rootfind = SciMLBase.LeftRootFind,
    reinitializealg = SciMLBase.CheckInit()
)

# regenerate a new ODESystem from the old @mtkmodel-generated ODESystem, but with the new downcrossing event callback added
substitutions = Dict(eq.lhs => eq.rhs for eq in observed(odesys_ref))
substituted_eqs = [substitute(eq, substitutions) for eq in equations(odesys_ref)] # make sure the external forcing time series is represented in the right place
@named odesys_new = ODESystem(substituted_eqs, get_iv(odesys_ref), unknowns(odesys_ref), parameters(odesys_ref), observed=observed(odesys_ref),  continuous_events=[cb_down])

odesys = structural_simplify(odesys_new);