Hello,
In a nutshell what I do is the following :
1- First I monitor the evolution of one output through a continuous callback and when it hits a limit value I saturate the derivative output by multiplying it by a boolean set to false .
2- Through a discrete callback that becomes active aftetwards, I compute(via interpolation) at each integration step the expected value at t+dt without saturation. If that value is below limit I deactivate the saturation.
My problem is that the value extrapolated when the saturation is deactivated does not evolve while it should. So I want to know if the integrator’s extrapolation is affected by the modification of the parameters it holds during a callback. Or what am I missing ?
There goes a reproducible example :
using DifferentialEquations
using Plots
#params
H = 0.0371
ΔPe = 0.05
saturated = [false]
observer = [false]
p = (H, ΔPe, saturated)
#DAE definition
function system!(du, u, p, t)
H, ΔPe, saturated = p
du[1] = -5.0u[1] - 1.0u[5]
du[2] = 0.0228532u[1] - 0.0263158u[2]
du[3] = 0.789474u[1] + 6.0u[2] -2u[3]
du[4] = (-0.263158u[1] -2.0u[2] +1.0u[3] - 1.0u[4] ) *!saturated[1]
du[5] = (u[4] -ΔPe)/(2*H)
end
#Continuous Callback
function condition(u, t, integrator)
u[4] - 0.06
end
function affect!(integrator)
observer[1] = true
integrator.p[3] .= true
integrator.u[4] = 0.06
end
cb1 = ContinuousCallback(condition, affect!, nothing);
#Discrete Callback
function condition_(u, t, integrator)
Bool(observer[1])
end
function affect_!(integrator)#Check if the value is increasing without any limitation; if yes, continue to saturate
integrator.p[3] .= false
s_t = integrator.sol(integrator.t)
s_tplus = integrator.sol(integrator.t+integrator.dt)#evaluate next value without saturation
@show s_t[4], s_tplus[4]
#
if s_tplus[4] < 0.06 #if under limit, deactivate saturation
observer[1] = false
elseif s_tplus[4] >= 0.06#if still above limit, maintain saturation
integrator.p[3] .= true
end
end
cb2 = DiscreteCallback(condition_, affect_!)
cb = CallbackSet(cb1, cb2);
#DAE Solving
M = [1.0 0 0 0 0
0 1.0 0 0 0
0 0 1.0 0 0
0 0 0 0 0
0 0 0 0 1.0]
u0 = [0.0, 0.0, 0.0, 0.0, 0.0]
tspan = (0.0,30.0)
#with_logger(TerminalLogger()) do
prob = ODEProblem(ODEFunction(system!, mass_matrix = M),u0,tspan, p)
sol = solve(prob,Rodas5(), reltol = 1e-8, abstol = 1e-8, callback = cb)
This is the output from @show :
(s_t[4], s_tplus[4]) = (0.06, 0.06)
(s_t[4], s_tplus[4]) = (0.06, 0.06)
(s_t[4], s_tplus[4]) = (0.06, 0.06)
(s_t[4], s_tplus[4]) = (0.06, 0.06)
