We need better tests on the behavior for this case, so please donate an example whether it’s determined to be working or not.
If you only have one VectorContinuousCallback then it’s not modified. If you have multiple in a CallbackSet, then it could be.
@ChrisRackauckas I’ve created #515 and #516 for some odd behavior I came across working this.
Here, I am going to show some other potential issues I came across but with a more involved example. I have a simple double integrator with a constant positive acceleration. It has an initial negative velocity. The system should get saturated at a lower limit. Once it is saturated, the dynamics for \dot{x} are set to zero. Additionally, b/c the system stops, \dot{x}=0.0 and the dynamics for \ddot{x} are set to zero as well. I am checking for three events
using Plots; gr()
using OrdinaryDiffEq
function eom!(du,u,params,t)
issaturated1, issaturated2 = params
du[1] = u[2] * !Bool(issaturated1)
du[2] = 3.0 * !any(Bool.([issaturated1, issaturated2]))
nothing
end
function condition!(out, u, t, integrator, limits)
out[3] = -1.0 #initialize to force 0 crossing on out[2] event detection
if integrator.p[1] == false #if withing saturation limits
out[1] = (u[1] - limits[1][2]) #up crossing of state 1 upper bound
out[2] = -(u[1] - limits[1][1]) #up crossing of state 2 lower bound
else
# turn saturation off and check for sign change of ẍ
u2 = deepcopy(u)
p2 = deepcopy(integrator.p)
p2[1:2] .= false
du = similar(u)
integrator.f(du,u2,p2,t)
out[3] = du[2]
end
@show out, t
end
function affect!(integrator, event_index)
@show event_index
if event_index ∈ [1,2]
# Set velocity to zero
for u_ in full_cache(integrator)
u_[2] = 0.0
end
# saturate state 1 and unsaturate state 2 b/c it goes to 0.0
integrator.p[1] = true
integrator.p[2] = false
end
end
begin
u0 = [0.0,-10.0]
tspan = (0.0,8.0)
# Set saturation limits, compute if IC in limits, and set params accordingly
saturationlimits = [[-10.0,Inf],[-Inf,Inf]]
issaturated = .!([lim[1] <= u_ <= lim[2] for (u_,lim) ∈ zip(u0, saturationlimits)])
params = [issaturated...]
prob_sat = ODEProblem(eom!,u0,tspan,params)
cb = VectorContinuousCallback((out, u, t, integrator) -> condition!(out, u, t, integrator,saturationlimits),
affect!,3, affect_neg! = nothing, save_positions=(true,true))
sol_nocb = solve(prob_sat,Tsit5()) # solve w/o callback
sol_cb = solve(prob_sat,Tsit5(), callback = cb) # solve w/ callback
end
begin
hline(saturationlimits[1],label="saturation limit",color=:black, ls=:dashdotdot)
hline!(saturationlimits[2],label="saturation limit",color=:black, ls=:dashdotdot)
plot!(sol_nocb,leg=false, ls=:dash)
plot!(sol_cb)
end
This produces the following error on sol_cb = solve(prob_sat,Tsit5(), callback = cb)
BoundsError: attempt to access 3-element view(::Array{Float64,1}, 1:3) with eltype Float64 at index [-1]
__subarray_throw_boundserror(::Type{SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}}, ::Array{Float64,1}, ::Tuple{UnitRange{Int64}}, ::Int64, ::Int64, ::Tuple{Int64}) at subarray.jl:45
throw_boundserror at subarray.jl:43 [inlined]
checkbounds at abstractarray.jl:449 [inlined]
getindex(::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}, ::Int64) at subarray.jl:260
find_callback_time(::OrdinaryDiffEq.ODEIntegrator{Tsit5,true,Array{Float64,1},Float64,Array{Bool,1},Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Bool,1},ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},DiffEqBase.DEStats},ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{VectorContinuousCallback{getfield(Main, Symbol("##86#88")),typeof(affect!),Nothing,typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64,DiffEqBase.CallbackCache{Array{Float64,N} where N,Array{Float64,N} where N}}, ::VectorContinuousCallback{getfield(Main, Symbol("##86#88")),typeof(affect!),Nothing,typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}, ::Int64) at callbacks.jl:537
find_first_continuous_callback at callbacks.jl:221 [inlined]
handle_callbacks!(::OrdinaryDiffEq.ODEIntegrator{Tsit5,true,Array{Float64,1},Float64,Array{Bool,1},Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Bool,1},ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},DiffEqBase.DEStats},ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{VectorContinuousCallback{getfield(Main, Symbol("##86#88")),typeof(affect!),Nothing,typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64,DiffEqBase.CallbackCache{Array{Float64,N} where N,Array{Float64,N} where N}}) at integrator_utils.jl:247
_loopfooter!(::OrdinaryDiffEq.ODEIntegrator{Tsit5,true,Array{Float64,1},Float64,Array{Bool,1},Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Bool,1},ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,N...
In this case, out[2] crosses zero, and p[1] is set to true. So, for the next condition call it should check if the acceleration crosses zero, which it does by design. However, affect! is never called with event_idx=3.
If I @show out, t in condition! I get the following. Note how out[3] crosses zero.
...
(out, t) = ([-Inf, 6.27587e-12, -1.0], 1.2251482265554061)
(out, t) = ([-Inf, -0.0, -1.0], 1.2251482265544138)
(out, t) = ([-Inf, -0.0, -1.0], 1.2251482265544138)
event_index = 2
(out, t) = ([-Inf, -0.741555, 3.0], 1.2251482265544138)
(out, t) = ([-Inf, -0.0, 3.0], 8.0)
(out, t) = ([-Inf, -0.0, 3.0], 1.9779095347150344)
...
Furthermore, if I change affect! so integrator.p[1] = false, effectively putting the event logic in affect! for event_index = 1 or 2. The code doesn’t error. And produces the desired result.
function affect!(integrator, event_index)
@show event_index
if event_index == 1
# Set velocity to zero
for u_ in full_cache(integrator)
u_[2] = 0.0
end
# saturate state 1 and unsaturate state 2 b/c it goes to 0.0
integrator.p[1] = false
integrator.p[2] = false
end
end
So, I decided to rethink the problem and put special logic where this type of domino triggering of events is avoided. Here, the logic is in affect! for hitting a lower bound. This works.
using Plots; gr()
using OrdinaryDiffEq
function eom!(du,u,params,t)
issaturated1, issaturated2 = params
du[1] = u[2] * !Bool(issaturated1)
du[2] = 3.0 * !any(Bool.([issaturated1, issaturated2]))
nothing
end
function condition!(out, u, t, integrator, limits)
if integrator.p[1] == false #if withing saturation limits
out[1] = (u[1] - limits[1][2]) #up crossing of state 1 upper bound
out[2] = -(u[1] - limits[1][1]) #up crossing of state 2 lower bound
end
end
function affect!(integrator, event_index)
@show event_index
if event_index ∈ [1,2]
# Set velocity to zero and unsaturate state 2, assumes 0 ∈ bounds of state 2
for u_ in full_cache(integrator)
u_[2] = 0.0
end
integrator.p[2] = false
u2 = deepcopy(integrator.u)
p2 = deepcopy(integrator.p)
p2[1:2] .= false
du = similar(u2)
integrator.f(du,u2,p2,integrator)
if event_index == 2
integrator.p[1] = du[2] > 0.0 ? false : true
elseif event_index == 1
integrator.p[1] = du[2] < 0.0 ? false : true
end
end
end
begin
u0 = [0.0,-10.0]
tspan = (0.0,8.0)
# Set saturation limits, compute if IC in limits, and set params accordingly
saturationlimits = [[-10.0,Inf],[-Inf,Inf]]
issaturated = .!([lim[1] <= u_ <= lim[2] for (u_,lim) ∈ zip(u0, saturationlimits)])
params = [issaturated...]
prob_sat = ODEProblem(eom!,u0,tspan,params)
cb = VectorContinuousCallback((out, u, t, integrator) -> condition!(out, u, t, integrator,saturationlimits),
affect!,2, affect_neg! = nothing, save_positions=(true,true))
sol_nocb = solve(prob_sat,Tsit5())
sol_cb = solve(prob_sat,Tsit5(), callback = cb)
end
begin
hline(saturationlimits[1],label="saturation limit",color=:black, ls=:dashdotdot)
hline!(saturationlimits[2],label="saturation limit",color=:black, ls=:dashdotdot)
plot!(sol_nocb,leg=false, ls=:dash)
plot!(sol_cb)
end
Now, if I modify affect! to include logic for hitting an upper bound as well and change the saturation limits , i.e.
saturationlimits = [[-10.0,10.0],[-Inf,Inf]]
function affect!(integrator, event_index)
@show event_index
if event_index ∈ [1,2]
# Set velocity to zero and unsaturate state 2, assumes 0 ∈ bounds of state 2
for u_ in full_cache(integrator)
u_[2] = 0.0
end
integrator.p[2] = false
u2 = deepcopy(integrator.u)
p2 = deepcopy(integrator.p)
p2[1:2] .= false
du = similar(u2)
integrator.f(du,u2,p2,integrator)
if event_index == 2
integrator.p[1] = du[2] > 0.0 ? false : true
elseif event_index == 1
integrator.p[1] = du[2] < 0.0 ? false : true
end
end
end
Then things don’t go so well. I get the following error on sol_cb = solve(prob_sat,Tsit5(), callback = cb)
ArgumentError: The interval [a,b] is not a bracketing interval.
You need f(a) and f(b) to have different signs (f(a) * f(b) < 0).
Consider a different bracket or try fzero(f, c) with an initial guess c.
init_state(::Roots.AlefeldPotraShi, ::Roots.DerivativeFree{getfield(DiffEqBase, Symbol("##354#355")){OrdinaryDiffEq.ODEIntegrator{Tsit5,true,Array{Float64,1},Float64,Array{Bool,1},Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Bool,1},ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},DiffEqBase.DEStats},ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{VectorContinuousCallback{getfield(Main, Symbol("##38#40")),typeof(affect!),Nothing,typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64,DiffEqBase.CallbackCache{Array{Float64,N} where N,Array{Float64,N} where N}},VectorContinuousCallback{getfield(Main, Symbol("##38#40")),typeof(affect!),Nothing,typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}}}, ::Tuple{Float64,Float64}) at bracketing.jl:449
#find_zero#4(::Roots.NullTracks, ::Bool, ::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:atol,),Tuple{Float64}}}, ::Function, ::Function, ::Tuple{Float64,Float64}, ::Roots.AlefeldPotraShi, ::Nothing) at find_zero.jl:533
#find_zero at none:0 [inlined]
#find_zero at none:0 [inlined]
find_callback_time(::OrdinaryDiffEq.ODEIntegrator{Tsit5,true,Array{Float64,1},Float64,Array{Bool,1},Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Bool,1},ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},DiffEqBase.DEStats},ODEFunction{true,typeof(eom!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{VectorContinuousCallback{getfield(Main, Symbol("##38#40")),typeof(affect!),Nothing,typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64,DiffEqBase.CallbackCache{Array{Float64,N} where N,Array{Float64,N} where N}}, ::VectorContinuousCallback{getfield(Main, Symbol("##38#40")),typeof(affect!),Nothing,typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}, ::Int64) at callbacks.jl:504
find_first_continuous_callback at callbacks.jl:221 [inlined]
handle_callbacks!(::OrdinaryDiffEq.ODEIntegrator{Tsit5,true,Array{Float64,1},Float64,Array{Bool,1},Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float6...
Interestingly though, if I run the debugger and step all the way through, the code runs without error. Could this potentially be due to floating point difference with the interpreted code?
sol_cb = @Juno.run solve(prob_sat,Tsit5(), callback = cb)
Sorry for the drawn out example. 
I hope it make sense.
Any further progress or conclusions?
I’m interested in this feature but I cannot find a cool way of saturating integrator.
This issue was fixed a few years ago with the callback improvements in 2019. You’ll need a new thread with an example.