Integral saturation limits w/ DifferentialEquations.jl

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:

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

Roger that. Travel safely.

@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

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)
            push!(saturated, true)

# 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

function condition_saturate(u,t,integrator, idx, limit)
    u[idx] - limit

function affect_saturate!(integrator, idx)
    for u_ in full_cache(integrator)
        u_.saturated[idx] = true

function condition_unsaturate(u,t,integrator, idx)
    # check for sign change of derivate.
    u2 = deepcopy(u)
    u2.saturated .= false
    du = similar(u)
    return du[idx]

function affect_unsaturate!(integrator, idx)
    for u_ in full_cache(integrator)
        u_.saturated[idx] = false

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)

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)

1 Like

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.

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.

Fixed in . We’ll get that into a release soon. Thanks for the example.

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.

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.

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?

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

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

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?

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

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?

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

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.

Because of the OP’s post I added a specific protection that prevents this from occurring:

DrPapa, do you have an example using these 2 callback functions cb1, cb2 using the “ison” boolean variable. I am trying to do the same thing as you- place integral saturation limits on states in for my model described by a system of ODEs. Thanks.

I do… somewhere. I will post once I can find it. Might be faster to make a new example :flushed:

1 Like

@Matthew_Overlin I have not forgotten about you. I found my example and played with it some more and found that it was really brittle and that the “ison” approach has some issues. Primarily the issue is with the non-zero value to return

Unfortunately, I don’t think this is true, the sign of the value does matter. For

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

The sign of the non-zero value needs to be consistent with func1(u) otherwise when ison[1] switches, an artificial zero-crossing will be detected. For Example, if you have ison ? -0.2 : 1, a zero-crossing will trigger when ison flips

I’m going to spend some more time looking at it to see if I can come up with a better solution.

1 Like

yeah, I guess that’s a caveat to my statement: it doesn’t matter as long as you get the sign right.

When this thread was originally posted VectorContinuousCallback wasn’t an option. After playing with it some, I think the need for using a Ref(ison) goes away. @ChrisRackauckas is out in condition!(out,u,t,integrator) ever modified externally to this function? If not, I think simply not updating the value in condition! will give the effect of switching the callback off and this sign issue goes away since you are effectively returning the last value computed.

I think it is working but I may have hit a bug in VectorContinuousCallback when you have multiple events in the same step or when an event immediately triggers another. In the later case the second is never detected even though the sign in out changes. It very well may be an issue on my end. I will try to put a MWE together today to test it.