Reachability using Taylor model flowpipes for nonlinear ODEs

[mforets]

Papers:

• Taylor Model Flowpipe Construction
  for Non-linear Hybrid Systems. Chen, Abraham, Sankaranarayanan
• Implementation of Taylor Models in CORA 2018. Althoff, Grebenyuk, Kochdumper.
• Using Taylor Models in the Reachability Analysis
  of Non-linear Hybrid Systems. Chen, Abraham, Sankaranarayanan.
• Rigorous integration of flows using Taylor Models. Makino, Berz.

Theses:

• Reachability Analysis of Non-Linear Hybrid Systems using Taylor Models. Xin Chen's PhD Thesis., in particular chapter 3.
• Mioara Joldes, Rigorous Polynomial Approximations and Applications
• Kyoko Makino, Rigorous Analysis of Nonlinear Motion in Particle Accelerators

Implementations:

• Flow *
• CORA

Julia packages:

• Taylor models are implemented in JuliaIntervals/TaylorModels.jl
• Taylor series are implemented in JuliaDiff/TaylorSeries.jl.

[roccaa]

The thesis of Chen! You should also put the link to Flow* code which implements this thesis.

[schillic]

Taylor models were recently implemented in CORA with a detailed description in this ARCH18 paper.

[mforets]

Thanks for the link. In Julia, Taylor models are implemented in JuliaIntervals/TaylorModels.jl and there's the related package JuliaDiff/TaylorSeries.jl.

[lbenet]

These are some further interesting links

• Mioara Joldes, Rigorous Polynomial Approximations and Applications

• Kyoko Makino, Rigorous Analysis of Nonlinear Motion in Particle Accelerators

I highly recommend Joldes' thesis.

[mforets]

In JuliaReach/Reachability.jl#537 we added the TMJets continuous post. The flowpipe is represented as a collection of hyperrectangles, and the algorithm is defined in validated_integ, at the moment in the validatedODEs branch. Thanks @lbenet 🌮 🌮 🎉 !

Eventually, it may be more convenient to just write the reachability algorithm in this repo (for extensions, eg. to do flowpipe-invariant intersections in the loop).

Examples of use are in this notebook

• https://nbviewer.jupyter.org/github/mforets/escritoire/blob/master/reachability/BenchsHybrid.ipynb?flush_cache=true

• https://github.com/mforets/escritoire/blob/master/reachability/BenchsNLN.ipynb and the folder ReachabilityBenchmarks/models/ARCH/NLN.

---

I think that we can close this issue, but before i would like to store the research literature links somewhere visible, maybe in the Reachability docs.

[lbenet]

Sorry for the late comment on this... Can you comment about the time for each computation?

[mforets]

If the modes are of the form x'(t) = Ax(t) + u(t), i would expect linear continuous posts to be more efficient than nonlinear ones for the same accuracy, since the linear structure can be exploited, particularly for high-dimensional systems. But the bouncing ball is "middle-ground" for different approaches, since it is only two-dimensional.

If you check again the BenchsHybrid notebook you'll see updated plots and times. At the moment the three algorithms are in the same order ~0.5-1.0sec for a time horizon of 12sec, which gives 3 jumps. Note also that some time is spent on the discrete post, which is the same function for all three. Another benchmark would be to compare the behavior without jumps.

---

Incidentally for the bouncing ball model with @taylorize i had to write

@taylorize function bball!(t, x, dx)
    dx[1] = x[2]
    dx[2] = -1.0 + 0.0001*x[1] # here 0.0 * x[1] doesn't work, see below
    return dx
end

i cannot pass dx[1] = -1.0 + 0.0 * x[1], dx[2] = -one(x[1]) or dx[2] = -1.0 + zero(x[1]) because this gives a δt in validated_step! that is too large (the whole time horizon T) such that solve only returns one set. I checked that validated_integ returns a single set. @lbenet Should i open an issue? In which package?

[mforets]

Another benchmark would be to compare the behavior without jumps.

OTOH the circle model, also 2D, see BenchsNLN notebook, takes ~10ms for BFFPSV and ~600ms for TMJets. It would be interesting to profile the bouncing ball to see how the time distribute between discretizaiton, continuous post and discrete post.

[lbenet]

Should i open an issue? In which package?

Please do so. I think it should be opened in TaylorIntegration.jl since @taylorize lives there.

[lbenet]

i cannot pass dx[1] = -1.0 + 0.0 * x[1], dx[2] = -one(x[1]) or dx[2] = -1.0 + zero(x[1]) because this gives a δt in validated_step! that is too large (the whole time horizon T) such that solve only returns one set. I checked that validated_integ returns a single set.

I've checking what is happening here. I think the following explains the issue: The equations of motion correspond to constant acceleration, so they are solved (exactly!) with a polynomial of degree 2 in t. If you use orderT larger than 3 you will see the problem. The problem arises because the time step is determined by the last two coefficients of the local solution, which for orderT larger than 3 are both identically zero; TaylorIntegration.stepsize! will then return a stepsize which is Inf, which is then reduced to tmax (time horizon) by TaylorIntegration.taylorstep!. In other words, the polynomial of order 2 in t is the exact solution, for all t!

We haven't come with a way to automatically consider previous coefficients, but I am thinking about this now.

@mforets Could you check that things behave "properly" if you set orderT=2 or orderT=3? Could you also open an issue at TaylorIntegration/jl ?

[mforets]

Thanks for having a close look into it. Actually with orderT equal to 2 or 3 it doesn't work: validated_integ returns a list of sets but they are all practically the same as the initial set. Also tried varying orderQ and abs_tol.

Could you also open an issue at TaylorIntegration/jl ?

See PerezHz/TaylorIntegration.jl#66

[lbenet]

I'll check it in more detail...