#718 - Add Flowpipe type #719

schillic · 2019-12-12T17:37:05Z

Closes #718.
Closes #707.

add changes
test changes (mainly projection and plotting, particularly for the operators whose project methods got removed) in practice

These changes turned out to be more complicated than I thought because we do have several places where we manually extract information from the vectors of reach sets (which are now Flowpipes) and we have a certain plotting feature that required more work.

Notable changes:

new type Flowpipe
ReachSolution now stores a Vector of Flowpipes.
- new helper function n_sets(::ReachSolution) that counts the number of sets
Continuous-post operators now return a Flowpipe object instead of a Vector{<:AbstractReachSet}.
I removed the additional project methods in some continuous-post operators (they were not callable anymore with the change of the interface).
tube⋂inv! does not receive the Vector of reach sets anymore but rather computes a new one; consequently I removed the ! in the name.
A bigger change was needed to account for plotting certain indices of a ReachSolution only. The code should still give in the previous result, but it is obviously more complicated now that the sets distribute over different flowpipes.

schillic · 2019-12-14T15:26:07Z

This is what I tested with. In master this crashes in various places. So this PR is already an improvement.

using Reachability, MathematicalSystems, HybridSystems, LinearAlgebra
using TaylorIntegration
using TaylorModels: @taylorize
using Reachability: solve
using LazySets: HalfSpace

A = [0.0509836  0.168159  0.95246   0.33644
     0.42377    0.67972   0.129232  0.126662
     0.518654   0.981313  0.489854  0.588326
     0.38318    0.616014  0.518412  0.778765];
X0 = BallInf(ones(4), 0.1);
B = [0.866688  0.771231
     0.065935  0.349839
     0.109643  0.904222
     0.292474  0.227857];
U = Interval(0.99, 1.01) × Interval(0.99, 1.01);

s1 = solve(IVP(LCS(A), X0), :T=>0.1);
s2 = solve(IVP(LCS(A), X0),
          Options(:T=>0.1, :project_reachset=>true),
          op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4]));
s3 = solve(IVP(LCS(A), X0),
          Options(:T=>0.1, :project_reachset=>false),
          op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4]));
s4 = project(s3);

s5 = solve(IVP(CLCCS(A, B, nothing, U), X0),
           Options(:T=>0.2, :project_reachset=>true), op=GLGM06());

@taylorize function vanderPol!(dx, x, params, t)
    local μ = 1.0
    dx[1] = x[2]
    dx[2] = (μ * x[2]) * (1 - x[1]^2) - x[1]
    return dx
end;
𝑆 = BlackBoxContinuousSystem(vanderPol!, 2);
X0 = Hyperrectangle(low=[1.25, 2.35], high=[1.55, 2.45]);
𝑃 = InitialValueProblem(𝑆, X0);
s6 = solve(𝑃, Options(:T=>7.0, :mode=>"reach"),
           op=TMJets(:abs_tol=>1e-10, :orderT=>10, :orderQ=>2));
s7 = solve(𝑃, Options(:T=>7.0, :mode=>"reach", :project_reachset=>true),
           op=TMJets(:abs_tol=>1e-10, :orderT=>10, :orderQ=>2));

function bouncing_ball()
    # automaton structure
    automaton = LightAutomaton(1)

    # mode 1
    A = [0.0 1.0; 0.0 0.0]
    B = reshape([0.0, -1.0], (2, 1))
    U = Singleton([1.0])
    inv = HalfSpace([-1.0, 0.0], 0.0) # x >= 0
    m1 = ConstrainedLinearControlContinuousSystem(A, B, inv, ConstantInput(U))

    # modes
    modes = [m1]

    # transition from m1 to m1 (self-loop)
    add_transition!(automaton, 1, 1, 1)
    A = [1.0 0.0; 0.0 -0.75]
    guard = HPolyhedron([HalfSpace([0.0, 1.0], 0.0),   # v ≤ 0
                         HalfSpace([-1.0, 0.0], 0.0),  # x ≥ 0
                         HalfSpace([1.0, 0.0], 0.0)])  # x ≤ 0
    t1 = ConstrainedLinearMap(A, guard)

    # transition annotations
    resetmaps = [t1]

    # switching
    switchings = [AutonomousSwitching()]

    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)

    # initial condition in mode 1
    X0 = Hyperrectangle(low=[10.0, 0.0], high=[10.2, 0.0])
    initial_condition = [(1, X0)]

    system = InitialValueProblem(ℋ, initial_condition)

    options = Options(:mode=>"reach", :T=>5.0,
                      :plot_vars=>[1, 2], :max_jumps=>1, :verbosity=>1)

    return (system, options)
end

system, options = bouncing_ball();
sols = [];
for opC in [BFFPSV18(:δ=>0.1), GLGM06(:δ=>0.1)]
    sol = solve(system, options, opC, LazyDiscretePost(:check_invariant_intersection => true));
    sol_proj = project(sol);
    push!(sols, sol_proj);
end

using Plots
plot(s2)
plot(s4)
plot(s5)
plot(s5, indices=[1,3,5,7,9,20])
plot(s6)
plot(s7)
plot(sols[1])
plot(sols[2])

mforets · 2019-12-14T15:46:17Z

This is what I tested with. In master this crashes in various places. So this PR is already an improvement.

Sounds good!

Does this PR also solve some or all of #707 ?

Also, as the PR introduces breaking changes i would like to tag a new minor version of this package before merging this PR, for future reference.

schillic · 2019-12-14T21:12:52Z

Does this PR also solve some or all of #707 ?

It resolves the crashes.

using Reachability, MathematicalSystems
A = rand(4, 4);
S = LinearContinuousSystem(A);
X0 = rand(Hyperrectangle, dim=4);
M = rand(1, 4);
P = InitialValueProblem(S, X0);
opts = Options(:T=>1.0, :project_reachset=>true, :projection_matrix=>M);
opts_algo = Options(:δ=>0.1, :vars=>[1, 2]);
solve(P, opts, op=BFFPSV18(opts_algo));

#707 uses inputs that are not supported (M = rand(2, 4)), but I would say that then an error is expected.

mforets · 2019-12-14T21:22:45Z

Good. So we can as well close #707 with this PR, what do you think?

mforets · 2019-12-14T21:26:12Z

This is what I tested with. In master this crashes in various places. So this PR is already an improvement.

Can you link to a notebook with the plots?

mforets · 2019-12-14T21:29:25Z

I'm curious because bouncing ball was not working well last time i checked. In this notebook i tried Reachability#master with the bouncing ball example and the overapproximation error was very bad. If you scroll down to Out [85] there is a picture that shows GLGM06 using zonotopes and with a smaller overapproximation error of the intersection.

schillic · 2019-12-15T08:07:04Z

Can you link to a notebook with the plots?

using Reachability, MathematicalSystems, HybridSystems, LinearAlgebra, Plots

function bouncing_ball()
    # automaton structure
    automaton = LightAutomaton(1)

    # mode 1
    A = [0.0 1.0; 0.0 0.0]
    B = reshape([0.0, -1.0], (2, 1))
    U = Singleton([1.0])
    inv = HalfSpace([-1.0, 0.0], 0.0) # x >= 0
    m1 = ConstrainedLinearControlContinuousSystem(A, B, inv, ConstantInput(U))

    # modes
    modes = [m1]

    # transition from m1 to m1 (self-loop)
    add_transition!(automaton, 1, 1, 1)
    A = [1.0 0.0; 0.0 -0.75]
    guard = HPolyhedron([HalfSpace([0.0, 1.0], 0.0),   # v ≤ 0
                         HalfSpace([-1.0, 0.0], 0.0),  # x ≥ 0
                         HalfSpace([1.0, 0.0], 0.0)])  # x ≤ 0
    t1 = ConstrainedLinearMap(A, guard)

    # transition annotations
    resetmaps = [t1]

    # switching
    switchings = [AutonomousSwitching()]

    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)

    # initial condition in mode 1
    X0 = Hyperrectangle(low=[10.0, 0.0], high=[10.2, 0.0])
    initial_condition = [(1, X0)]

    system = InitialValueProblem(ℋ, initial_condition)

    options = Options(:mode=>"reach", :T=>20.0,
                      :plot_vars=>[1, 2], :max_jumps=>3, :verbosity=>1)

    return (system, options)
end

system, options = bouncing_ball();
sols = [];
for opC in [BFFPSV18(:δ=>0.01), GLGM06(:δ=>0.01)]
    sol = solve(system, options, opC, LazyDiscretePost(:check_invariant_intersection => true));
    sol_proj = project(sol);
    push!(sols, sol_proj);
end

plot(sols[1])
plot(sols[2])

mforets · 2019-12-15T11:01:27Z

The plots look good. I can't tell if the second plot has zonotopes, is it the case?

schillic · 2019-12-15T11:02:37Z

I can't tell if the second plot has zonotopes, is it the case?

The algorithm is GLGM06, so I think those are zonotopes.

mforets · 2019-12-15T11:14:32Z

Ok, i was asking because that didn't work on master, see #716.

mforets · 2019-12-15T11:15:12Z

If the PR solves that problem that issue can be closed as well.

schillic · 2019-12-15T11:21:49Z

Ok, i was asking because that didn't work on master, see #716.

I see. No, this PR did not address that issue.

julia> flowpipe_GLGM06 = sols[2].flowpipes[1];

julia> typeof(flowpipe_GLGM06.reachsets[2])
ReachSet{Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}}

mforets · 2019-12-15T11:24:58Z

Ok, thanks.

schillic added 6 commits December 12, 2019 11:09

add Flowpipe type

9235c29

use Flowpipe type

ab9bbaf

rewrite project helpers

490c4a5

remove project helpers

2b55c0a

refactor interface of solve and discrete post

00201ca

add test that projects

c1dea38

schillic force-pushed the schillic/718 branch from 302cee7 to c1dea38 Compare December 14, 2019 14:44

schillic changed the title ~~WIP #718 - Add Flowpipe type~~ #718 - Add Flowpipe type Dec 14, 2019

schillic requested a review from mforets December 14, 2019 15:25

mforets approved these changes Dec 14, 2019

View reviewed changes

schillic mentioned this pull request Dec 15, 2019

Errors with projection #707

Closed

schillic merged commit 573612f into master Dec 15, 2019

schillic deleted the schillic/718 branch December 15, 2019 11:09

schillic mentioned this pull request Jan 23, 2020

WIP Remove type parameters from ReachSolution #694

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

#718 - Add Flowpipe type #719

#718 - Add Flowpipe type #719

schillic commented Dec 12, 2019 •

edited

Loading

schillic commented Dec 14, 2019

mforets commented Dec 14, 2019

schillic commented Dec 14, 2019

mforets commented Dec 14, 2019

mforets commented Dec 14, 2019

mforets commented Dec 14, 2019

schillic commented Dec 15, 2019

mforets commented Dec 15, 2019

schillic commented Dec 15, 2019

mforets commented Dec 15, 2019 •

edited

Loading

mforets commented Dec 15, 2019

schillic commented Dec 15, 2019

mforets commented Dec 15, 2019

#718 - Add Flowpipe type #719

#718 - Add Flowpipe type #719

Conversation

schillic commented Dec 12, 2019 • edited Loading

schillic commented Dec 14, 2019

mforets commented Dec 14, 2019

schillic commented Dec 14, 2019

mforets commented Dec 14, 2019

mforets commented Dec 14, 2019

mforets commented Dec 14, 2019

schillic commented Dec 15, 2019

mforets commented Dec 15, 2019

schillic commented Dec 15, 2019

mforets commented Dec 15, 2019 • edited Loading

mforets commented Dec 15, 2019

schillic commented Dec 15, 2019

mforets commented Dec 15, 2019

schillic commented Dec 12, 2019 •

edited

Loading

mforets commented Dec 15, 2019 •

edited

Loading