@ChrisRackauckas I have a “working” example below. I put working in quotes b/c if the time span is increased (e.g. 20.0s) I get an error: Double callback crossing floating pointer reducer errored. Report this issue. in find_callback_time at DiffEqBase/ZQVwI/src/callbacks.jl:264 I am on Julia 1.0 w/ DiffEqBase v5.5.1
Here is the results up to 18s.
Here is the result via Simulink
In this example, I make a SaturatedType{T} <: DEDataVector{T} type that keeps track of the saturation limits and a bool indicating whether the state is saturated. I then end up making 3 ContinuousCallbacks for each state (upper limit, lower limit, and call back on derivative of the unsaturated dynamics to check when to turn saturation off).
As I was writing this I ultimately wanted to be able to turn callbacks off via callbacks. E.g. when the upper limit callback triggers for state 1 I would like to be able to switch that callback off and then switch on the “unsaturate” callback for that state. This would essentially enable a simple finite state machine via callbacks.
Do you have any thoughts on a better way to approach this problem?
using Plots; gr()
using DifferentialEquations
##
mutable struct SaturatedType{T} <: DEDataVector{T}
x::Array{T,1} # state
limits::Array{Array{T,1}} # Upper/Lower Saturation Limits
saturated::Array{Bool,1} # Indicator if saturated == true
end
function SaturatedType(x,limits)
@assert length(x) == length(limits)
# Loop through each state and determine value for saturated
saturated = Bool[]
for (val, lim) ∈ zip(x,limits)
if lim[1] <= val <= lim[2]
push!(saturated, false)
else
push!(saturated, true)
end
end
SaturatedType(x,limits,saturated)
end
# Spring Mass Damper w/ saturation
function smd_saturated!(du,u,params,t)
m,k,c = params
du[1] = u[2]
du[2] = -(k*u[1]+c*u[2])/m + 20*sin(t)
du .= du .* .!u.saturated #saturate
nothing
end
function condition_saturate(u,t,integrator, idx, limit)
u[idx] - limit
end
function affect_saturate!(integrator, idx)
for u_ in full_cache(integrator)
u_.saturated[idx] = true
end
end
function condition_unsaturate(u,t,integrator, idx)
# check for sign change of derivate.
u2 = deepcopy(u)
u2.saturated .= false
du = similar(u)
integrator.f(du,u2,integrator.p,t)
return du[idx]
end
function affect_unsaturate!(integrator, idx)
for u_ in full_cache(integrator)
u_.saturated[idx] = false
end
end
function make_callbacks(limits)
# Loop through each state and create callbacks for upper/lower saturation
# and for "unsaturating" the integrator
n = length(limits)
cbl = ContinuousCallback[]
for ii ∈ 1:n
affect_sat = (integrator) -> affect_saturate!(integrator, ii)
affect_unsat = (integrator) -> affect_unsaturate!(integrator, ii)
condition_lower = (u,t,integrator) -> condition_saturate(u,t,integrator,ii, limits[ii][1])
condition_upper = (u,t,integrator) -> condition_saturate(u,t,integrator,ii, limits[ii][2])
condition_unsat = (u,t,integrator) -> condition_unsaturate(u,t,integrator,ii)
cb_lower = ContinuousCallback(condition_lower,nothing, affect_neg! = affect_sat)#, idxs=ii)
cb_upper = ContinuousCallback(condition_upper,affect_sat, affect_neg! = nothing)#, idxs=ii)
cb_unsat = ContinuousCallback(condition_unsat, affect_unsat)
push!(cbl, cb_lower, cb_upper, cb_unsat)
end
cbl
end
save_positions = (true,true)
u0 = [0.0,0.0]
limits = [[-5.0,5.0],[-10.0,10.0]]
u0_sat = SaturatedType(u0,limits) #make "saturated" IC type
tspan = (0.0,18.0)
params = [1.0, 1.0, 0.0]
cbs = CallbackSet(make_callbacks(limits)...)
prob_sat = ODEProblem(smd_saturated!,u0_sat,tspan,params)
sol = solve(prob_sat,Tsit5(),callback = cbs)
##
plot(sol)
