JSoC 2019: Improved flowpipe/guard intersections for hybrid reachability using Taylor model

Hi there!

I'm Aadesh and was accepted for JSoC2019. My project is to add new flowpipe/guard intersection algorithms and improve the way this is handled to have more accuracy and scalability in reachability analysis applications. My mentors are Marcelo Forets, Luis Benet, and David Sanders. This is my first blog post as I started my project, I'd like to share both technical and general aspects throughout the summer.

What is Reachability?

In a dynamical system, given an initial condition, one can find out the motion of the system by integrating the ordinary differential equation (ODE). The goal of flowpipe approximation is to build an enclosure of the ODE given a set of initial conditions and uncertain inputs or uncertain parameters, thus obtaining the sets of state reachable by the dynamical system. All trajectories are guaranteed to be contained in the computed flowpipe. One such set-based method is the Taylor models approach that we consider in this project, for the case of hybrid dynamical systems. An example of hybrid system is given in the next paragraph.

Example of bouncing ball with air friction

My first contribution in JuliaReach was to add a new hybrid system test case (see ReachabilityBenchmarks#134), the nonlinear bouncing ball from (1).

The bouncing ball can be viewed as a hybrid system equipped with two state variables and which represent the vertical distance from the ball to the ground and the velocity of the ball, respectively (see the picture below). The velocity evolves both continuously when the ball is moving in the air and discretely when the ball is bouncing up. The model has two discrete states and discrete jumps.

The model has two guards:

: 0 0 : −0.8· denotes that if the value of is zero and has a negative value then the jump may be executed and value of is updated,
0 implies that when velocity becomes zero the ball starts to fall down.

using Pkg
Pkg.add("MathematicalSystems")
Pkg.add("LazySets")
Pkg.add("LinearAlgebra")
Pkg.add("HybridSystems")
Pkg.add("Plots")
Pkg.clone("https://github.com/JuliaReach/Reachability.jl.git")
Pkg.add("TaylorIntegration")

54.5s

Julia

In our implementation of this example, bball_up and bball_down are two ODEs representing the motion of the ball as it goes up and down respectively. The function bouncing_ball captures the whole process, and returns the associated initial value problem for this hybrid system. We will consider uncertain initial conditions in position and in speed:

the position is [4.9, 5.1]
the velocity is [−0.2, 0]

As in the above figure we have two modes, invariant, transitions and Guard each defining the flow of the process. the detailed description of all the terms can be found here.

233.1s

Julia

using MathematicalSystems, Reachability, LinearAlgebra, HybridSystems
using LazySets
using Reachability: solve
using Plots
using TaylorIntegration


@taylorize function bball_up!(t, x, dx) 
    dx[1] = x[2]
    dx[2] = -9.8 - 0.1*(x[1])^2
    return dx
end

@taylorize function bball_down!(t, x, dx)
    dx[1] = x[2]
    dx[2] = -9.8 + 0.1*(x[1])^2
    return dx
end

function bouncing_ball()

    automaton = LightAutomaton(2) # two nodes

    inv_up = HPolyhedron([HalfSpace([-1.0, 0.0], 0.0),  # x >= 0
                          HalfSpace([0.0, -1.0], 0.0)]) # v >= 0

    inv_down = HPolyhedron([HalfSpace([-1.0, 0.0], 0.0),# x >= 0
                         HalfSpace([ 0.0, 1.0], 0.0)])  # v <= 0

    m1 = ConstrainedBlackBoxContinuousSystem(bball_up!, 2, inv_up)

    m2 = ConstrainedBlackBoxContinuousSystem(bball_down!,2, inv_down)

    modes = [m1,m2]  #modes

    add_transition!(automaton, 2, 1, 1)    #alpha transition
    add_transition!(automaton, 1, 2, 2)    #beta transition



    guard_alpha = HPolyhedron([HalfSpace([1.0, 0.0], 0.0),   #x>=0
                              HalfSpace([-1.0, 0.0], 0.0)])  #x<=0
    guard_beta =  HPolyhedron([HalfSpace([0.0, 1.0], 0.0),   #v<=0
                               HalfSpace([0.0,-1.0], 0.0)])  #v>=0


    A = [1.0 0.0; 0.0 -0.75]
    t1 = ConstrainedLinearMap(A, guard_alpha)
    t2 = ConstrainedLinearMap(A, guard_beta)

    #resetmaps
    resetmaps = [t1,t2]

    # switching
    switching = AutonomousSwitching()
    switchings = fill(switching, 2)

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

    #initial condition in mode one
    X0= Hyperrectangle(low=[ 4.9,-0.2], high=[ 5.1,0.0])

    initial_condition = [(2,X0)]

    system = InitialValueProblem(ℋ, initial_condition)

    options = Options(:mode=>"reach", :T=>6.0, :plot_vars=>[1, 2],                                                        :project_reachset=>false)

    return (system, options)
end


problem, options = bouncing_ball();

options = Options(:mode=>"reach", :T=>3.0, :plot_vars=>[1, 2],                                :project_reachset=>false, :verbosity => "info")

@time sol_TMJets = solve(problem, options, TMJets(:orderT=>5, :orderQ=>2,      :abs_tol=>1e-10),LazyDiscretePost(:check_invariant_intersection=>true));


plot(sol_TMJets, use_subindices=false, aspectratio=0.5, alpha=.5,                                                            color=:lightblue)

233.1s

Julia

Let me clear some things here, Reachability works well in this case but as the number of dimension and complex nature increases, we get big box instead of the smooth plot one can check it by making small changes into the ODE's used for the Bouncing Ball model. There are two main related difficulties: the wrapping effect and how intersections between the flowpipe and the guard are handled.

Wrapping Effect

Algorithms for the approximation of reachable sets for nonlinear ODEs suffer because of "wrapping effect". To illustrate the wrapping effect, suppose that we start with an initial set given by a box and apply successive rotations. We will see that if we use box approximations after each operation, the size of the set grows.Consider the blue square in the picture below and apply a counter-clockwise rotation of = /4 obtaining the pink square.

using Pkg
Pkg.add("LazySets")
Pkg.add("Plots")
using LazySets, Plots

# box (ball in the infinity norm)
B = BallInf([2.5,1.0], 0.5)
plot(B)
# rotation of θ rad
R(θ) = [cos(θ) -sin(θ); sin(θ) cos(θ)]
plot!(R(pi/4)*B, ratio=1)

21.6s

Julia

Now suppose that we wrap the rotated square (pink) using an box (green) the same type as the original square (blue).

using LazySets.Approximations: box_approximation
const □ = box_approximation
plot!(□(R(pi/4)*B), ratio=1)
plot!(R(pi/4) * □(R(pi/4)*B), ratio=1)

2.1s

Julia

We see that subsequent approximations of the rotated sets result in loose bounds. Using other set representations can be helpful to mitigate this problem.

plot!(□(R(pi/4) * □(R(pi/4)*B)), ratio=1)

1.9s

Julia

Groundwork

As discussed with my mentors I decided to start with some groundwork. It includes conversion of a Taylor model into a Zonotope (or a Polytope) so that we can have some strong guess how it's intersection with a guard will look like, whereas benchmarking those with some standard polynomial tell about the accuracy and the speed. Another approach is to split the domain so that tighter bounds can be obtained.

Over-approximate TaylorModel with Zonotope

I would like to start with converting a TaylorModel with linear terms only, we began to experiment with it since the conversion to a zonotope is exact. A Taylor model (TM) can be thought of as the tuple (P, I) where p is a polynomial and I is an interval remainder. For the example we consider a vector of two TMs constructed using the polynomials (p1, p2) and the interval remainder I.

using Pkg
Pkg.add("TaylorModels")
Pkg.add("LazySets")
Pkg.add("LinearAlgebra")
Pkg.add("Plots")
Pkg.add("IntervalArithmetic")

25.9s

Julia

In this method, we linearise the polynomial which simply means to remove the non-linear term from the polynomial part of the TM.

using Plots
using LazySets:Zonotope,BallInf
using TaylorModels
using IntervalArithmetic:Interval 
using LinearAlgebra:Diagonal
δ = 0.5; I = Interval(-δ, δ)
x₀ = Interval(0.0)
D = Interval(-3.0, 1.0)
p1 = Taylor1([2.0, 1.0], 2)
p2 = Taylor1([0.9, 3.0], 2)
TM = [TaylorModel1(p1, I, x₀, D), TaylorModel1(p2, I, x₀, D)]

154.3s

Julia

2-element Array{TaylorModel1{Float64,Float64},1}: 2.0 + 1.0 t + [-0.5, 0.5] 0.9 + 3.0 t + [-0.5, 0.5]

The first step is to normalize the polynomial (), in normalization the reformulate

() into () and [−1, 1] such that the ranges of the and remain same.

l_TM = length(TM)
Q = [normalize_taylor(TM[i].pol, D, true) for i = 1:l_TM]

2.4s

Julia

2-element Array{Taylor1{Float64},1}: 1.0 + 2.0 t + 𝒪(t³) - 2.1 + 6.0 t + 𝒪(t³)

Now collect the center and generators of the Zonotope () over approximating TM. Centre of the is the constant part of the and generator is the coefficient of which is th term in the

rem = [TM[i].rem for i = 1:l_TM]
α = [mid(rem[i]) for i = 1:l_TM]
rem_gen = [abs(rem[i].hi - α[i]) for i = 1:l_TM]#generator
c = [Q[i].coeffs[1] + α[i] for i = 1:l_TM] #centre 
gen = [Q[i].coeffs[2]  for i = 1:l_TM] #generator
z = Zonotope(c, hcat(gen, Diagonal(rem_gen)));plot(z)

43.1s

Julia

Box over (Pink) approximation is just the product of Range of each polynomial as a set and the Remainder interval is added to each of them.

e = evaluate(TM[1].pol, D) +  TM[1].rem
r = evaluate(TM[2].pol, D) +  TM[2].rem
K = e × r
plot!(K)

3.3s

Julia

Finally, look for the accuracy by plotting the actual TM, as you can see in the following result the green plot is the TM evaluation and the zonotope(Blue) is tightly over-approximates the TM, so this method works very well.

using LazySets:Singleton,⊕
D = Interval(-3.0, 1.0)
B = BallInf([0.0, 0.0], 0.5) 
fig = plot(z)
fig = plot!(fig, K)
for ti in D.lo:0.1:D.hi
    x = 2.0 + 1.0*ti
    y = 0.9 + 3.0*ti

    plot!(fig, Singleton([x, y]) ⊕ B, legend=:none, aspectratio=.5,                                                               color="green")
end
fig

7.4s

Julia

Time for TM with non-linear polynomial. Step of creating 'TM' is the same besides the polynomial (P1, P2) we use have non-linear terms.

using Pkg
Pkg.add("TaylorModels")
Pkg.add("LazySets")
Pkg.add("LinearAlgebra")
Pkg.add("Plots")
Pkg.add("IntervalArithmetic")

72.0s

Julia

using Plots
using LazySets:Zonotope,BallInf
using TaylorModels
using IntervalArithmetic:Interval 
using LinearAlgebra:Diagonal
p1 = Taylor1([2.0, 1.0, 1.0], 2)
p2 = Taylor1([0.9, 3.0], 2)
δ = 0.5; I = Interval(-δ, δ)
x₀ = Interval(0.0)
D = Interval(-3.0, 1.0)
TM = [TaylorModel1(p1, I, x₀, D),
      TaylorModel1(p2, I, x₀, D)]

160.8s

Julia

2-element Array{TaylorModel1{Float64,Float64},1}: 2.0 + 1.0 t + 1.0 t² + [-0.5, 0.5] 0.9 + 3.0 t + [-0.5, 0.5]

In this method, we linearise the polynomial which simply means to remove the non-linear term.

l_TM = length(TM)
pol_lin = [constant_term(TM[i].pol) + linear_polynomial(TM[i].pol) for i =                                                                1:l_TM]
rem_r = [TM[i].rem for i = 1:l_TM]

2.0s

Julia

2-element Array{Interval{Float64},1}: [-0.5, 0.5] [-0.5, 0.5]

Then just take Box over-approximate of non-linear part.

pol_nonlin = [TM[i].pol - pol_lin[i] for i = 1:l_TM]
rem_nonlin = [evaluate(pol_nonlin[i], D) for i = 1:l_TM] 
rem_nonlin = rem_r .+ rem_nonlin

1.1s

Julia

2-element Array{Interval{Float64},1}: [-3.5, 9.5] [-0.5, 0.5]

Q = [normalize_taylor(TM[i].pol, D, true) for i = 1:l_TM]

0.6s

Julia

2-element Array{Taylor1{Float64},1}: 2.0 - 2.0 t + 4.0 t² + 𝒪(t³) - 2.1 + 6.0 t + 𝒪(t³)

The plot below is the zonotope (Blue) over-approximation of the TM,

α = [mid(rem_nonlin[i]) for i = 1:l_TM]
rem_gen = [abs(rem_nonlin[i].hi - α[i]) for i = 1:l_TM]
c = [Q[i].coeffs[1] + α[i] for i = 1:l_TM] #center 
gen = [Q[i].coeffs[2]  for i = 1:l_TM] #generator
Z = Zonotope(c, hcat(gen, Diagonal(rem_gen)));plot(Z)

1.5s

Julia

In a similar way, we take the box over-approximation (Pink).

e = evaluate(TM[1].pol, D) + TM[1].rem
r = evaluate(TM[2].pol, D) + TM[2].rem
T = e×r
plot!(T)

1.3s

Julia

Actual TM (green) as you can see in the following plot and the zonotope (Blue) is poor over approximates of TM, so this doesn't work. again it might help in other ways like if one sees the over approximate as the set of Zonotope then using set operations can get the better results.

using LazySets:Singleton,⊕
fig = plot(Z)
fig = plot!(fig, T)
B = BallInf([0.0, 0.0], 0.5)
for ti in D.lo:0.1:D.hi 
     x = 2 + 1.0*ti + 1.0*ti^2
     y = 0.9 + 3.0*ti
     plot!(fig, Singleton([x, y]) ⊕ B, legend=:none, aspectratio=.5,                                                              color="green")
end 
fig

1.2s

Julia

What Next?

We are working the over-approximation of a Taylor model using support functions (polytopic overapproximation) which is an optimization problem at its heart. There are a few ways to approach it, it is under discussion now. I have also implemented a "Branch and Bound" algorithm and I would like to compare it with the other existing methods and software. Special thanks to Christian Schilling JuliaReach member for help. With these, I'd like to conclude.

sources