JuliaReach / ReachabilityAnalysis.jl

Computing reachable states of dynamical systems in Julia
https://juliareach.github.io/ReachabilityAnalysis.jl/
MIT License
189 stars 17 forks source link

Reachability using the parallelotope method #93

Open mforets opened 5 years ago

mforets commented 5 years ago

See:

dpsanders commented 4 years ago

@mforets Do you have an implementation of parallelotopes?

mforets commented 4 years ago

Hi David, A set type specific to parallelotopes is not available yet in LazySets, but we are working on it in LazySets#1632 -- PR which is kind of ready..

Depending on the application that you have in mind, you can use zonotopes, which are a generalization of parallelotopes of order greter than one.

As explained in the notes, there are different approaches for the parallelotope, either in constraint form (as in #1632) on in generator form. Again, depending on the use case, one representation may be preferable than the other.

dpsanders commented 4 years ago

OK great, thanks. I was asking in terms of implementing the algorithms in this and related papers by Goldsztejn.

mforets commented 4 years ago

Alright, I had a look to that paper again, and to get started we would need Theorem 3.1 in page 7, which gives a parallelotopic enclosure of the image set of a parallelotope through a (differentiable) function.

I have to think about it, but perhaps i can add it as an extension in RangeEnclosures.jl.

mforets commented 4 years ago

To mimic the notation in the paper, one can use a LazySets.AffineMap, which can be created with

Parallelotope(A, ๐ฎ, xฬƒ) = AffineMap(A, ๐ฎ, xฬƒ)

or equivalently

Parallelotope(A, ๐ฎ, xฬƒ) = xฬƒ โŠ• A*๐ฎ

After looking with some detail Section 3.2, I think there are typos. Since xฬƒ is a vector, then ๐ฒ := ๐›š(xฬƒ) (mind the bold omega, for the interval extension of the discrete map omega) can only possibly be ฯ‰(xฬƒ), because by definition of interval extension, ฯ‰(x) = ๐›š(x) if x is a scalar. But then ๐ฒ - yฬƒ is identically zero, from Eq. (9). This doesn't make sense. In sum, shouldn't Theorem 3.1 define ๐ฒ := ๐›š(โ–กโŸจ๐ฑโŸฉ)?

mforets commented 4 years ago

With such considerations, I tried implementing Theorem 3.1 and the result can be found in this notebook. The final result is at the very end. I'm under the impression that i did something wrong, because the parallelotopic overapproximation is much bigger than the box overapproximation, which can be found in the plot in Out[31].

dpsanders commented 4 years ago

That's a really fantastic notebook, thanks! It illustrates very well a lot of the concepts that you've implemented so nicely in LazySets.jl.

I'll try to think about what's going on.

dpsanders commented 4 years ago

I agree that there is an error as you point out in the paper, and I think your interpretation is correct. I also agree that the result is a huge over-estimation. I think that's because the original parallelotope is so big -- I believe these estimates are supposed to be good for small initial objects?