# Integral saturation limits w/ DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143>\
**Category:** General Usage\
**Tags:** question, diffeq\
**Created:** [March 21, 2019, 1:45pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143 "2019-03-21T13:45:15Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [March 21, 2019, 1:45pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/1 "2019-03-21T13:45:16Z")

</div>

I need to simulate an ODE w/ saturation limits on the integrators to prevent integral wind-up, e.g. see [link](https://www.mathworks.com/help/simulink/slref/integrator.html#f7-890704.f7-1739821). 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 😉

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 22, 2019, 2:52pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/2 "2019-03-22T14:52:17Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [March 22, 2019, 3:11pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/3 "2019-03-22T15:11:23Z")

</div>

Roger that. Travel safely.

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [March 26, 2019, 4:50pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/4 "2019-03-26T16:50:07Z")

</div>

@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.

 ![11%20PM](https://global.discourse-cdn.com/julialang/original/3X/1/1/11cdc9b64b47b4c9e2cf6e23ea4d4ad60fd609d0.jpeg)

Here is the result via Simulink

 ![2019-03-26-11-07-57](https://global.discourse-cdn.com/julialang/original/3X/9/e/9ec8e34f9eff505fa12e97c3ad06a2ad12203045.png)

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?

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [March 26, 2019, 5:43pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/5 "2019-03-26T17:43:34Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 13, 2019, 11:35pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/6 "2019-04-13T23:35:31Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 14, 2019, 1:33am UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/7 "2019-04-14T01:33:31Z")

</div>

Fixed in [Safer epsilon in callback and handle callback switch off by ChrisRackauckas · Pull Request #207 · SciML/DiffEqBase.jl · GitHub](https://github.com/JuliaDiffEq/DiffEqBase.jl/pull/207) . We’ll get that into a release soon. Thanks for the example.

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [April 15, 2019, 2:02pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/8 "2019-04-15T14:02:37Z")

</div>

> [@ChrisRackauckas](#):
>
> I didn’t forget about this just had to get caught up.

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 15, 2019, 2:04pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/9 "2019-04-15T14:04:21Z")

</div>

> [@DrPapa](#):
>
> 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.

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [April 15, 2019, 2:25pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/10 "2019-04-15T14:25:32Z")

</div>

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

```julia
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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 16, 2019, 3:33am UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/11 "2019-04-16T03:33:17Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![korsbo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/korsbo/32/15942_2.png) [@korsbo](https://discourse.julialang.org/u/korsbo)\
**Post date:** [April 18, 2019, 4:31pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/12 "2019-04-18T16:31:47Z")

</div>

> [@ChrisRackauckas](#):
>
> 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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 18, 2019, 4:35pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/13 "2019-04-18T16:35:46Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![korsbo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/korsbo/32/15942_2.png) [@korsbo](https://discourse.julialang.org/u/korsbo)\
**Post date:** [April 18, 2019, 4:55pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/14 "2019-04-18T16:55:27Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 18, 2019, 5:05pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/15 "2019-04-18T17:05:13Z")

</div>

> [@korsbo](#):
>
> 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: [https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/callbacks.jl#L205](https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/callbacks.jl#L205)

---

<div class="post-metadata">

**Author:** ![Matthew\_Overlin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/matthew_overlin/32/10533_2.png) [@Matthew\_Overlin](https://discourse.julialang.org/u/Matthew_Overlin)\
**Post date:** [October 1, 2019, 7:37pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/16 "2019-10-01T19:37:48Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [October 2, 2019, 3:44pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/17 "2019-10-02T15:44:58Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [October 10, 2019, 4:48pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/18 "2019-10-10T16:48:28Z")

</div>

@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

> [@ChrisRackauckas](#):
>
> It doesn’t matter what the non-zero value is.

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

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

```

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 10, 2019, 9:35pm UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/19 "2019-10-10T21:35:49Z")

</div>

> [@DrPapa](#):
>
> 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

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

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [October 11, 2019, 11:40am UTC](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143/20 "2019-10-11T11:40:04Z")

</div>

> [@ChrisRackauckas](#):
>
> 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.

[Next page](https://discourse.julialang.org/t/integral-saturation-limits-w-differentialequations-jl/22143.md?page=2)