Integral saturation limits w/ DifferentialEquations.jl

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

  1. Does the state x come in contact with the upper bound, if so turn saturation on
  2. Does the state x come in contact with the lower bound, if so turn saturation on
  3. If x is saturated on the lower bound AND the acceleration becomes +'ve then reset the dynamics to the original EOM.
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)

sat

Sorry for the drawn out example. :grimacing: I hope it make sense.