ERROR: LoadError: Non bracketing interval passed in bisection method. Please report the error in DiffEqBase. #599
Description
I have a simple julia script which implements an example with saturation limits on ODE states. I am trying to see one state saturate on the lower limit. I am able to see a state saturate at the upper limit. I received an error of "Non bracketing interval passed in bisection method" and I was told to report this error. To reproduce the issue, I'm seeing, run the julia code below. The example I'm running here is minorly changed from the example from DrPapa here: https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/22
import PyPlot;
#using Plots; gr()
using DifferentialEquations;
# Taken from DrPaPa @:
# https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/22
# VectorContinuousCallback documentation:
# https://diffeq.sciml.ai/stable/features/callback_functions/
# The actual system of ODEs. This is a simple system of ODEs here.
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
# This condition function is continuously checked throughout
# the whole time-domain simulation.
function condition!(out, u, t, integrator, limits)
if integrator.p[1] == false #if state u1 is within 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 1 lower bound
end
end
# The "affect!" function below seems to only be called
# when either out[1] or out[2] is zero in the condition! function above.
# That is, when either the state u1 has reached a lower or upper saturation
# limit, this "affect!" function is called. printing "event_index" will show
# you if out[1] was zero or out[2] was zero.
function affect!(integrator, event_index)
println("\ninside affect! function..");
@show event_index
if event_index ∈ [1,2] # 1 of the 2 sat. limits was met.
# Set velocity to zero and unsaturate state 2, assumes 0 ∈ bounds of state 2
println("integrator: ", integrator);
println("full_cache(integrator): ", full_cache(integrator));
println("typeof(full_cache(integrator)): ", typeof(full_cache(integrator)));
println("integrator.p: ", integrator.p);
println("integrator.u: ", integrator.u);
println("integrator.du: ", integrator.du);
println("integrator.t: ", integrator.t);
#println("size(full_cache(integrator)): ", size(full_cache(integrator)));
# Set velocity to zero and unsaturate state 2, assumes 0 ∈ bounds of state 2
for u_ in full_cache(integrator)
#println("u_: ", u_);
u_[2] = 0.0;
#println("u_: ", u_);
end
#println("u_: ", u_);
println("full_cache(integrator): ", full_cache(integrator));
println("integrator.p: ", integrator.p);
#integrator.p[2] = false;
println("integrator.p: ", integrator.p);
u2 = deepcopy(integrator.u);
println("u2: ", u2);
p2 = deepcopy(integrator.p);
p2[1:2] .= false
du = similar(u2)
integrator.f(du,u2,p2,integrator); # this evaluate the integrator ???
if event_index == 2 # u1 reaches lower limit.
integrator.p[1] = du[2] > 0.0 ? false : true
elseif event_index == 1 # u1 reaches upper limit.
integrator.p[1] = du[2] < 0.0 ? false : true
end
println("integrator.p: ", integrator.p);
end
end
function affect2!(integrator, event_index)
println("\ninside affect! function..");
@show event_index
if event_index ∈ [1,2] # 1 of the 2 sat. limits was met.
# Set velocity to zero and unsaturate state 2, assumes 0 ∈ bounds of state 2
println("integrator: ", integrator);
println("full_cache(integrator): ", full_cache(integrator));
println("typeof(full_cache(integrator)): ", typeof(full_cache(integrator)));
println("integrator.p: ", integrator.p);
println("integrator.u: ", integrator.u);
println("integrator.du: ", integrator.du);
println("integrator.t: ", integrator.t);
integrator.p[2] = false;
println("integrator.p: ", integrator.p);
u2 = deepcopy(integrator.u)
u2[2] = 0.0;
p2 = deepcopy(integrator.p)
p2[1:2] .= false
du = similar(u2)
integrator.f(du,u2,p2,integrator); # this evaluate the integrator ???
if event_index == 2 # u1 reaches lower limit.
integrator.p[1] = du[2] > 0.0 ? false : true
elseif event_index == 1 # u1 reaches upper limit.
integrator.p[1] = du[2] < 0.0 ? false : true
end
#integrator.p[1] = true;
println("integrator.p: ", integrator.p);
end
end
# ----------------------------------------
# The actual simulation is below.
# ----------------------------------------
time_to_sim = 8.0;
u0 = [0.0,-10.0];
tspan = (0.0,time_to_sim);
data_times = collect(0.0:1.0e-3:time_to_sim);
# Set saturation limits, compute if IC in limits, and set params accordingly
saturationlimits = [[-10.0,10.0],[-Inf,Inf]];
issaturated = .!([lim[1] <= u_ <= lim[2] for (u_,lim) ∈ zip(u0, saturationlimits)])
params = [issaturated...];
#cb = VectorContinuousCallback((out, u, t, integrator) -> condition!(out, u, t, integrator, saturationlimits),
# affect2!,2, affect_neg! = nothing, save_positions=(true,true));
cb = VectorContinuousCallback((out, u, t, integrator) -> condition!(out, u, t, integrator, saturationlimits),
affect2!,2, save_positions=(true,true));
prob_sat = ODEProblem(eom!,u0,tspan,params);
sol_nocb = solve(prob_sat,Tsit5(), saveat=data_times);
sol_cb = solve(prob_sat,Tsit5(), callback = cb, saveat=data_times);
sol_nocb_t = sol_nocb.t;
matrix_nocb = Array(sol_nocb);
u1_nocb = matrix_nocb[1,:];
u2_nocb = matrix_nocb[2,:];
sol_cb_t = sol_cb.t;
matrix_cb = Array(sol_cb);
u1_cb = matrix_cb[1,:];
u2_cb = matrix_cb[2,:];
# Plots
PyPlot.close("all");
PyPlot.figure();
PyPlot.plot(sol_cb_t, saturationlimits[1][1] .* ones(length(sol_cb_t)), "-.");
PyPlot.plot(sol_cb_t, saturationlimits[1][2] .* ones(length(sol_cb_t)), "-.");
PyPlot.plot(sol_nocb_t, u1_nocb, "--");
PyPlot.plot(sol_nocb_t, u2_nocb, "--");
PyPlot.plot(sol_cb_t, u1_cb);
PyPlot.plot(sol_cb_t, u2_cb);
PyPlot.xlabel("Time [sec.]");
PyPlot.legend(["u1 low limit", "u1 up limit", "u1_nocb", "u2_nocb", "u1_cb", "u2_cb"], loc="best");