Recording entry and exit time of state space sets with DifferentialEquations.jl

Specific Domains Modelling & Simulations
1

I am trying to estimate entry and exit times in sets of the state space, as functionality for DynamicalSystems.jl. I have a working implementation for discrete systems but for continuous it turns out it is much harder. Or hopefully not, perhaps you guys have a better idea.

The figure bellow describes my situation:

image
image708×566 29 KB

I have an ODE that I evolve in time and it gives the continuous black line (in theory) and the blue points in reality. I define some state u_0, and I need to record the times that the trajectory enters a ball of size ε_1 centered in u_0, and also the times that the trajectory exits this ball. The figure shows several u_0 purely for visualization, in my application only 1 u_0 will exist.

I also want to do this computation for several ball sizes ε_1, ε_2, …, at once, for performance reasons: if the trajectory cannot enter ball of size ε_1 (e.g it is in state m+1 in the figure), then it absolutely cannot enter ε_2 and thus the check can be skipped.

Now, I am really not sure how to go about this. Initially I wanted to use the callback functionality of diff eq, to minimize the amount of work I need to do. But I am not sure how to extract t_1 and t_2 (and save each in a different container, because the first is entry times, the second is exit times).

EDIT HERE:

The callbacks provide directly the interpolated time

and therefore I can directly store this time in some array. Therefore a simple way to do it (but probably not most performant) is to use ContinuousCallback.

2

You could use linear interpolation between each u in the trajectory and check to see if the minimum distance between the interpolated line and the center of the ball is equal to the radius of the ball.

1 Like
3

This was my original idea. I am a bit afraid because I don’t know whether the integrators can do interpolation like so:

image

i.e. that there is a loop that could be in principle be interpolated, but falls out of the linear connection.

4

The computed trajectory is discrete though, and if the true trajectory has additional features that the integrator doesn’t capture, you just need a better integrator.

1 Like
5

The loop in your picture would be populated by points that you could interpolate between

6

A different approach is to formulate the problem as a hybrid system with modes at the ball’s positions and identity reset maps. Then solve it using set-based integration and count the amount of time spent on each mode. This method is by construction mathematically rigorous and exhaustive (ie. covers all trajectories of the continuous system).

1 Like
7

Although I got to admit that I understood almost nothing of what you said, I should point out that I am interested in the return time to the ε-ball, not the time spent inside it.

8

Treating as a hybrid system is using ContinuousCallbacks like we discussed on the #diffeq-bridged channel. You’d just have an upcrossing affect! that fills one array and a downcrossing affect! that fills a different one. That should be a lot cheaper than the interval-based method, but the interval-based method has more of a guarantee (though in this kind of system I wouldn’t see it missing points).

Integrators already handle this: this is exactly the reason for the interp_points argument in callbacks, which defaults to 10 meaning it checks at 10 intermediate points of the trajectory for a zero crossing and thus would catch this (and most people never change this argument).

9

Alright, I am now very pleased to say that I have a working version for a single radius, that goes like this:

using DynamicalSystemsBase
using OrdinaryDiffEq: Tsit5
using DynamicalSystemsBase.DiffEqBase: ContinuousCallback, ODEProblem, solve
using Distances

εdistance(u, u0, ε::Real) = euclidean(u, u0) - ε

function transit_time_statistics(ds::ContinuousDynamicalSystem, u0, εs, T;
        alg = Tsit5(), diffeq...
    )
    eT = eltype(ds.t0)
    check_εs_sorting(εs, length(u0))
    exits = [eT[] for _ in 1:length(εs)]
    entries = [eT[] for _ in 1:length(εs)]

    # Make the magic callback:
    crossing(u, t, integ) = ChaosTools.εdistance(u, u0, εs[1])
    negative_affect!(integ) = push!(entries[1], integ.t)
    positive_affect!(integ) = push!(exits[1], integ.t)
    cb = ContinuousCallback(crossing, positive_affect!, negative_affect!;
        save_positions = (false, false)
    )

    prob = ODEProblem(ds, (eT(0), eT(T)); u0 = u0)
    sol = solve(prob, alg;
        callback=cb, save_everystep = false, dense = false,
        save_start=false, save_end = false, diffeq...
    )
    return exits, entries
end

This works so far without problem. You will notice that εs is a vector of ball radii, so I need to extend this framework to use a CallbackSet. So I change my callback generating code to:

    callbacks = ContinuousCallback[]
    for i in eachindex(εs)
        crossing(u, t, integ) = ChaosTools.εdistance(u, u0, εs[i])
        negative_affect!(integ) = push!(entries[i], integ.t)
        positive_affect!(integ) = push!(exits[i], integ.t)
        cb = ContinuousCallback(crossing, positive_affect!, negative_affect!;
            save_positions = (false, false)
        )
        push!(callbacks, cb)
    end
    cb = CallbackSet(callbacks...)

which, sometimes, works excellently. There are edge cases that I haven’t been able to identify yet where the integrator and/or interpolation seems to get stuck when exiting the ε-ball and or entering the inner most ε-ball and the callback positive_affect! gets called non-stop.

For example if I have εs = sort!([1.0, 0.1, 0.01]; rev=true), I can get resulting in output arrays like (fully reproducible code at the end):

julia> exits[2]
885033-element Array{Float64,1}:
     0.009521130925491618
   105.1285060629573
   227.79608608298875
   456.0488702006206
   561.549220491485
   954.0609666798003
  1610.4113223742547
  1645.4402235743025
  1680.468316963271
  1949.7593549776773
  2658.4544078395106
     ⋮

 27017.53815060816
 27017.53815060816
 27017.53815060816
 27017.53815060816
 27017.53815060816
 27017.53815060816

while the corresponding entries[2], which, by the definition of continuity of the trajectory, can have at most 1 more entry that exits[2], is:

julia> entries[2]
80-element Array{Float64,1}:
   105.10947207509459
   227.7840886145673
   456.03680102588226
   561.5319922039274
     ⋮
 23077.10816055895
 24143.055169575884
 24992.380689658865
 25009.884596962012
 25062.475812273526
 26355.939325594318
 26508.239300245827
 26783.405955508315
 27017.519471225707

So I guess something happens when exiting one of the inner balls. @ChrisRackauckas do you have any suggestions? In the example above I attempted to integrate for total of T = 100000.0 so the algorithm shouldn’t stop at 27017.519471225707 ( I guess it hit the max number of steps and thus terminated)

I can also say that I’ve noticed that when this stuck happens, it was the first time that my trajectory crossed the innermost ball of radius εs[3], which has the smallest radius.

Reproducible code (but unfortunatley, needs fixodprob branch of DynamicalSystemsBase)

using DynamicalSystemsBase
using OrdinaryDiffEq: Tsit5
using DynamicalSystemsBase.DiffEqBase: ContinuousCallback, ODEProblem, solve
using Distances

function transit_time_statistics(ds::ContinuousDynamicalSystem, u0, εs, T;
        alg = Tsit5(), diffeq...
    )
    eT = eltype(ds.t0)
    check_εs_sorting(εs, length(u0))
    exits = [eT[] for _ in 1:length(εs)]
    entries = [eT[] for _ in 1:length(εs)]

    # Make the magic callbacks:
    callbacks = ContinuousCallback[]
    for i in eachindex(εs)
        crossing(u, t, integ) = ChaosTools.εdistance(u, u0, εs[i])
        negative_affect!(integ) = push!(entries[i], integ.t)
        positive_affect!(integ) = push!(exits[i], integ.t)
        cb = ContinuousCallback(crossing, positive_affect!, negative_affect!;
            save_positions = (false, false)
        )
        push!(callbacks, cb)
    end
    cb = CallbackSet(callbacks...)

    prob = ODEProblem(ds, (eT(0), eT(T)); u0 = u0)
    sol = solve(prob, alg;
        callback=cb, save_everystep = false, dense = false,
        save_start=false, save_end = false, diffeq...
    )
    return exits, entries
end

ro = Systems.roessler(ones(3))
tr = trajectory(ro, 5000; Ttr = 100)
u0 = trajectory(ro, 11; Ttr = 100)[end] # return center
εs = sort!([1.0, 0.1, 0.01]; rev=true)
T = 100000.0 # maximum time
exits, entries = transit_time_statistics(ro, u0, εs, T)
1 Like
10

It should throw a warning. You might want to increase maxiters. You can add that to the ODEProblem you generate so it’s higher by default, since yes with that high of a T you might accidentally hit it with the defaults.

11

That is the problem: the integrator gets stuck at that point and never progresses. But I am not sure why. Increasing maxiters doesn’t help here, I clearly shouldn’t be getting the value 27017.53815060816 more than once.

12

Oh it might be hitting double callback… what DiffEqBase version do you have?

(@kanav99 is fixing some of these issues)

13 
(@v1.5) pkg> st DiffEqBase
Status `C:\Users\datse\.julia\environments\v1.5\Project.toml`
  [2b5f629d] DiffEqBase v6.44.3

In the mean time, could you please point me to the point of the source where the interpolation and calling of Roots.jl happens? I want to see if it is possible for me to somehow incorporate this information in my callbacks, so that I let DiffEq handle the first callback (outer most ε ball), but do the inner ones manually. Perhaps this will be both more stable and more performant?

14

You’ll want to update. This issue is likely from us using eps(typeof(t)) instead of eps(t) which was fixed in Callback bug fixes by pengwyn · Pull Request #574 · SciML/DiffEqBase.jl · GitHub . That should fix it.

The source for this is all in https://github.com/SciML/DiffEqBase.jl/blob/master/src/callbacks.jl#L650-L673 . Note that @kanav99 is currently writing a new rootfinder to fix the performance issues of Roots.jl

15

Hm, I did update to .47 but the problem remained the same. I’ll clean it up to a pure DiffEqBase MWE and open up an issue.