SciML
ContinuousEvent affect broken in mtkcompile?
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[])
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. 😃
