# ContinuousCallback for early termination of ODE fails to trigger

**URL:** <https://discourse.julialang.org/t/continuouscallback-for-early-termination-of-ode-fails-to-trigger/116910>\
**Category:** New to Julia\
**Tags:** diffeq, ode\
**Created:** [July 11, 2024, 1:29pm UTC](https://discourse.julialang.org/t/continuouscallback-for-early-termination-of-ode-fails-to-trigger/116910 "2024-07-11T13:29:40Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![apazo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apazo/32/216122_2.png) [@apazo](https://discourse.julialang.org/u/apazo)\
**Post date:** [July 11, 2024, 1:29pm UTC](https://discourse.julialang.org/t/continuouscallback-for-early-termination-of-ode-fails-to-trigger/116910/1 "2024-07-11T13:29:40Z")

</div>

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:

```julia
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.

```julia
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.

---

<div class="post-metadata">

**Author:** ![eldee](https://avatars.discourse-cdn.com/v4/letter/e/b5a626/32.png) [@eldee](https://discourse.julialang.org/u/eldee)\
**Post date:** [July 25, 2024, 11:02am UTC](https://discourse.julialang.org/t/continuouscallback-for-early-termination-of-ode-fails-to-trigger/116910/2 "2024-07-25T11:02:42Z")

</div>

TL;DR: Use `DiscreteCallback` instead of `ContinuousCallback`.  
  
  
Full answer:  
  
Hi apazo,

Welcome to the Julia community!

It seems to me that the problem is that you’re using a `ContinuousCallback` in the first place. Reading through the [documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/), a `ContinuousCallback` is triggered when a continuous function has a zero-crossing. As we are always working with discrete timesteps, this requires some interpolation. The benefit is that you can then work back to find the ‘precise’ time of the zero-crossing.  
If you want to use `ContinuousCallback`, your `condition_stop` should therefore return a float, not a boolean (though `Bool`s seem to get converted to `Float64` (`0.` or `1.`), I guess). Potentially you might be able to rewrite `condition_stop` to use something like

```julia
if ...
    return abs(integrator.p.saved_maxima[end][1] - integrator.p.saved_maxima[end - 1][1]) - integrator.p.tolerance # - instead of <
end

```

which changes from positive to negative when we should stop. But I couldn’t get this to work. In any case, your stopping criterion seems fundamentally discrete to me, as you’re explicitly relying on a discrete sampling `saved_maxima`. (As a continuous alternative, maybe you can check that the derivative is sufficiently close to 0 instead?)

If I simply keep your code for `condition_stop`, but change

```julia
cb_stop = ContinuousCallback(condition_stop, end_computation!)

```

to

```julia
cb_stop = DiscreteCallback(condition_stop, end_computation!)

```

I get as output

```julia
Stability target (ϵ ≤ 1.0e-5) has been reached at time 10.456329814160023 s, timestep #19

```

i.e. the stopping criterion does now trigger.

For reference, with your second (implementation of the) stopping criterion as a `DiscreteCallback`, I get

```julia
Stability target (ϵ ≤ 1.0e-5) has been reached at time 9.0 s, timestep #18

```

I imagine this occurs one timestep earlier, as in the iteration where you get a maximum,

```julia
(t in integrator.p.ON) & (integrator.p.I == integrator.p.I_max)

```

is already `true`, whereas in the

```julia
(integrator.sol.u[end - 1][1] > integrator.sol.u[end][1]) & (integrator.sol.u[end - 1][1] > integrator.sol.u[end - 3][1])

```

approach, this will only become `true` in the next iteration, when the maximum is found in `integrator.sol.u[end - 1][1]`.

(By the way, if I understand correctly, `t in integrator.p.ON` and `integrator.p.I == integrator.p.I_max` are equivalent, so you only need one of them? Also, you’d want to use `&&` (logical and) instead of `&` (bitwise and).)

What is unclear to me, though, is why `ContinuousCallBack` still works for the second implementation, but not for the first. If you look at the `return`ed signals, in the first implementation you get `false, false, ..., false, true, true, ..., true`, whereas in the second you get `false, false, ..., false, true, false, false`. So I guess the double ‘zero-crossing’ in the second case causes the callback to actually trigger. You can sort of simulate this in the first implementation by using

```julia
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 rand() < 0.5 # or rand() - 0.5
end 

```

(though you really shouldn’t). Perhaps someone who knows how `ContinuousCallBack` works under the hood can explain why `false, true` (`0.`, `1.`) does not count as a ‘zero-crossing’ (nor `true, false`, or `-0.5, 0.5` for that matter) while `false, true, false` does.

In any case, if you just use `DiscreteCallBack` you don’t need to worry about any of this 🙂.