Hi all,
I’m new to Julia, and I’m trying to simulate a simple ODE problem with some callbacks. My model consists of one single equation with two terms, one positive for growth, and one negative for decrease. However, the term for decrease depends on an external parameter I, which acts as a cyclic ON – OFF switch, the duration of an on-off cycle is defined by τ, and the fraction of the cycle where I is on is determined by duty. To represent the switch between on and off I’ve calculated all the timepoints where there would be a change in I, store those in 3 vectors (discontinuities, ON and OFF) and used the vectors to set up two discrete callbacks, where I save the states right before and after the callbacks.
My issue arises with the last callback I want to implement. As I follows a cyclic ON – OFF pattern, my state u ends up having a sawtooth pattern, with multiple local maxima and minima, but in the long term, its overall behavior is a decrease from its initial value of 1.0 to a certain stability point. I want to set up a continuous callback that stops the integration when the solution is close enough to that stable point, and for that, I set up a continuous callback that saves the local maxima and then compares the values of the last two saved maxima, and if they are within a certain tolerance, stops the integration.
I’ve run the script in the debugger from VSCode and verified that the continuous callback is indeed properly saving the values of the local maxima, and it does compare the value of the last maximum with the previous one, but for some reason it fails to trigger the callback. I leave you here a minimal working example of my code:
using DifferentialEquations, Plots, Sundials
I = 200.0
I_max = 200.0
I_min = 0.0
tolerance = 0.00001
discontinuities::Vector{Float64} = []
ON::Vector{Float64} = []
OFF::Vector{Float64} = []
saved_maxima::Vector{Tuple{Float64, Float64}} = []
timestep = 0.0
t_end = 45000
τ = 1.5*3
duty = 2/3
while timestep < t_end
push!(discontinuities, timestep + duty*τ, timestep + τ)
push!(OFF, timestep + duty*τ)
push!(ON, timestep + τ)
global timestep += τ
end
p = (; I, I_max, I_min, tolerance, discontinuities, ON, OFF, saved_maxima)
function ODE!(du, u, p, t)
growth = 10.0^(-4) * (1 - u[1])
decrease = 15.0^(-7) * p.I * u[1]
du[1] = growth - decrease
nothing
end
condition_turn_ON(u, t, integrator) = t in integrator.p.ON
function turn_ON!(integrator)
integrator.p = (; integrator.p..., I = integrator.p.I_max)
end
condition_turn_OFF(u, t, integrator) = t in integrator.p.OFF
function turn_OFF!(integrator)
integrator.p = (; integrator.p..., I = integrator.p.I_min)
end
function condition_stop(u, t, integrator)
# Check when a certain stability ϵ ( |u_i - u_i-1| < ϵ ) is reached
if length(integrator.sol.u) > 3
# Avoid out-of-bounds errors
if (integrator.sol.u[end - 1][1] > integrator.sol.u[end][1]) & (integrator.sol.u[end - 1][1] > integrator.sol.u[end - 3][1])
# Check whether a value is larger than the previous / next points. Maxima will be located at discontinuties, where data is saved twice, so [end - 2] = [end - 1]
if length(integrator.p.saved_maxima) == 0
# If it is the first maxima found, save it directly
push!(integrator.p.saved_maxima, (integrator.sol.u[end - 1][1], length(integrator.sol.u)))
elseif length(integrator.sol.u) != integrator.p.saved_maxima[end][2]
# If it is not the first maxima, check the length of integrator.sol to verify that it has been updated (it does not update every step like integrator.t)
push!(integrator.p.saved_maxima, (integrator.sol.u[end - 1][1], length(integrator.sol.u)))
end
end
end
if length(integrator.p.saved_maxima) >= 2
# Avoid out-of-bounds errors
if abs(integrator.p.saved_maxima[end][1] - integrator.p.saved_maxima[end - 1][1]) < integrator.p.tolerance
# Check if the absolute difference between two consecutive maxima is below the desired tolerance
return true
end
end
return false
end
function end_computation!(integrator)
println("Stability target (ϵ ≤ $(integrator.p.tolerance)) has been reached at time $(integrator.t) s, timestep #$(length(integrator.sol.t))")
terminate!(integrator)
end
cb_ON = DiscreteCallback(condition_turn_ON, turn_ON!, save_positions = (true, true))
cb_OFF = DiscreteCallback(condition_turn_OFF, turn_OFF!, save_positions = (true, true))
cb_stop = ContinuousCallback(condition_stop, end_computation!)
cbs_ODE = CallbackSet(cb_ON, cb_OFF, cb_stop)
ODE_model = ODEProblem(ODE!, [1.0], (0.0, t_end), p)
sol = solve(ODE_model, CVODE_BDF(), callback = cbs_ODE, tstops = discontinuities, progress = true)
NOTE: While writing this post I thought of a different way to define the continuous callback (I leave it below), this implementation saves the same maxima values and uses the same code to check if two maxima are close enough, but this one does trigger the callback successfully. Although it works, I`d still like to understand why is the first one failing.
saved_maxima::Vector{Float64} = []
function condition_stop(u, t, integrator)
if (t in integrator.p.ON) & (integrator.p.I == integrator.p.I_max) # maxima occur in the transtion from ON to OFF
push!(integrator.p.saved_maxima, integrator.sol.u[end - 1][1])
if length(integrator.p.saved_maxima) >= 2
if abs(integrator.p.saved_maxima[end] - integrator.p.saved_maxima[end - 1]) < integrator.p.tolerance
return true
end
end
end
return false
end
Thanks for your help, and best regards.
.