Integral saturation limits w/ DifferentialEquations.jl

question
#1

I need to simulate an ODE w/ saturation limits on the integrators to prevent integral wind-up, e.g. see link. My initial thought is to implement this via ContinuousCallback where the offending state’s dynamics are set to 0 during event detection. “Coming off” the limit would be handled by a separate callback that triggers when the saturation limit is met and the sign of the original dynamics changes appropriately.

Is this the right approach or does anybody have thoughts on a better way to tackle this problem? Calling @ChrisRackauckas :wink:

#2

I’m travelling so my systems are a bit off. Ping me in a day if I forget to respond.

#3

Roger that. Travel safely.

#4

@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)
1 Like
#5

Note that in the above example, the actual physics relating u[1] and u[2] are being ignored when applying the saturation. B/c it is a double integrator system u[2] should be reset to 0 anytime u[1] becomes saturated. It wouldn’t be hard to implement this manually in the affect functions, but it would be nice to automatically generate the callbacks.

#6

I didn’t forget about this just had to get caught up. Looking at it over the next few days and will get an answer. This issue is deep into floating point so I’m trying to find an adequate answer that is more accurate but doesn’t rely on too many rules about the user’s underlying code.

#7

Fixed in https://github.com/JuliaDiffEq/DiffEqBase.jl/pull/207 . We’ll get that into a release soon. Thanks for the example.

#8

No worries. I will give the fix a try and will let you know if I have any more issues. Thanks a lot.

Do you have any thoughts on this:

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.

Is this possible already and I am just missing something? If not, how feasible would it be, e.g. add/remove callbacks to a set during an event or turn callbacks in a set on/off during an event.

#9

I would use a Ref bool value to flip the callback. This has been requested enough though that we might want to add this to the interface.

#10

So, lets say I have 2 callbacks I need to switch between, are you suggesting the following:

function condition1(u,t,integrator, ison)
    ison[1] ? func1(u)  :    # some default non-zero value?
end

function condition2(u,t,integrator, ison)
    ison[2] ? func2(u)  :    # some default non-zero value?
end

function affect(integrator, ison)
    # update integrator
    ison .= .!ison
end

ison = [true, false]
cb1 = ContinuousCallback((u,t,integrator)->condition1(u,t,integrator, Ref(ison)), 
                                            integrator -> affect(integrator, Ref(ison))
cb2 = ContinuousCallback((u,t,integrator)->condition2(u,t,integrator, Ref(ison)), 
                                            integrator -> affect(integrator, Ref(ison))

if so, does it matter what is returned by the condition in the else branch, other then needing to be non-zero?

#11

Yes, that’s it. It doesn’t matter what the non-zero value is.

#12

Would you not get into trouble with ContinuousCallback if you all of a sudden switched sign in your condition function? Would it not try (and fail) to rootfind?

#13

I am assuming it’s a constant after it’s switched off

#14

I’m assuming so too, that is not what I meant.

I do, however, see my own mistake now. My worry was that if your condition function at one timestep returns a positive value and the next a negative value, then it would start trying to rootfind. But, in this particular example, the choice of return value will not cause a new root, since the previous step will, by design, have been a root itself.

So, if I’m correct, you could for example get into trouble if you used the same toggle for several different callbacks. Then the setting of your conditions return value to a constant might have caused it to initiate a rootfinding event where no actual root exists.

#15

Because of the OP’s post I added a specific protection that prevents this from occurring: https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/callbacks.jl#L205