ContinuousEvent affect broken in mtkcompile?

[hexaeder]

I am tyring to implement an anti-windup integrator using MTK and continuous events. This is the MWE:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as Dt
using OrdinaryDiffEqTsit5

@mtkmodel SaturatedIntegrator begin
    @parameters begin
        outMax=1 # Upper limi)t
        sat=0
    end
    @variables begin
        in(t), [description="Input signal"]
        out(t)=0, [description="Internal integrator state"]
        forcing(t)
    end
    @equations begin
        in ~ 10
        forcing ~ in - out
        Dt(out) ~ (1-sat) * forcing
    end
    @continuous_events begin
        [0 ~ out - outMax] => [sat ~ 1]
    end
end
@mtkcompile standalone = SaturatedIntegrator()
tspan = (0.0, 10.0)
prob = ODEProblem(standalone, Pair[], tspan)

The idea: as soon as out hits outMax, the parameter sat gets set to 1 which will disable the integrator. This MWE obviously doesn't include any unsaturation yet.
However, if i solve it, i get

julia> sol = solve(prob, Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 20-element Vector{Float64}:
  0.0
  9.999999999999999e-5
  0.0010999999999999998
  0.011099999999999997
  0.07603162341424605
  0.10536050712152709
  0.10536050712152709
  0.2533244579774494
  0.48384271375321186
  0.7693514972741362
  1.141672475399872
  1.5921004436502364
  2.1390680998900287
  2.7854045606631486
  3.5489527575506803
  4.445472246702398
  5.504531721512698
  6.765298284277327
  8.290366841448257
 10.0
u: 20-element Vector{Vector{Float64}}:
 [0.0]
 [0.0009999500016666248]
 [0.01099395221772342]
 [0.1103862230737223]
 [0.7321310206586767]
 [1.0]
 [1.0]
 [2.237840144532584]
 [3.8358984737102655]
 [5.366865322097515]
 [6.807152546164175]
 [7.9650182616803775]
 [8.822344153951112]
 [9.38293939989155]
 [9.712421399573932]
 [9.882641288156872]
 [9.95926428000569]
 [9.98841163887544]
 [9.99743504281877]
 [9.999515307769002]

where the double occurance of u=[1.0] implies, that the condition was met, however the affect didn't do anything.

I think the problem is, that the affect function was fundamentaly broken by mtkcompile, or at least I am confused by why the compiled model shows the entiere system equations as affect

julia> only(ModelingToolkit.get_continuous_events(standalone))
SymbolicContinuousCallback:
Conditions:
1-element Vector{Equation}:
 0 ~ -outMax + out(t)
Affect:
Affect system defined by equations:
Equation[in(t) ~ 10.0, forcing(t) ~ -out(t) + in(t)]
Negative-edge affect:
Affect system defined by equations:
Equation[in(t) ~ 10.0, forcing(t) ~ -out(t) + in(t)]

while before mtkcompile, the symbolic callback looks like i would've expected:

julia> @named nocompile = SaturatedIntegrator();
julia> only(ModelingToolkit.get_continuous_events(nocompile))
SymbolicContinuousCallback:
Conditions:
1-element Vector{Equation}:
 0 ~ -outMax + out(t)
Affect:
ModelingToolkit.SymbolicAffect(Equation[sat ~ 1], Equation[], Any[])
Negative-edge affect:
ModelingToolkit.SymbolicAffect(Equation[sat ~ 1], Equation[], Any[])

[zalanzsiboras]

Hi!

If you would like to change the sat parameter then you need to define it as a time dependent parameter: sat(t).

However, based on my testing this alone does NOT resolve the issue:

# Same with time dependent parameter
@mtkmodel SaturatedIntegrator begin
    @parameters begin
        outMax = 1 # Upper limit
        sat(t) = 0
    end
    @variables begin
        in(t), [description="Input signal"]
        out(t) = 0, [description="Internal integrator state"]
        forcing(t)
    end
    @equations begin
        in ~ 10
        forcing ~ in - out
        D(out) ~ (1 - sat) * forcing
    end
    @continuous_events begin
        [0 ~ out - outMax] => [sat ~ 1]
    end
end
@mtkcompile sys1 = SaturatedIntegrator()

prob1 = ODEProblem(sys1, [], (0.0, 10.0))
sol1 = solve(prob1)

It seems that the @continuous_events macro does not set sat as discrete parameter as it should.
Using SymbolicContinuousCallback() and building the system with ODESystem results in the correct solution:

function saturated_integrator_system()
@parameters outMax, sat(t)
@variables in_signal(t), out(t), forcing(t)
eqs = [
    in_signal ~ 10
    forcing ~ in_signal - out
    D(out) ~ (1 - sat) * forcing
]
event = [ModelingToolkit.SymbolicContinuousCallback([out ~ outMax] => [sat ~ 1], discrete_parameters = [sat])]

@mtkcompile sys = ODESystem(eqs, t, [in_signal, out, forcing], [outMax, sat]; continuous_events = event, defaults = [out => 0.0, outMax => 1.0, sat => 0])
end
sys2 = saturated_integrator_system()

prob2 = ODEProblem(sys2, Pair[], tspan)
sol2 = solve(prob2)

Unfortunately using SymbolicContinuousCallback() within @continuous_events is currently not possible, because it results in an error as I wrote in #3928 .

Of course it would be nice, if one of the maintainers would confirm that this is a bug. 😃

[ChrisRackauckas]

Yes this is a bug in the macro.