# 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

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")

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.

using MathematicalSystems, Reachability, LinearAlgebra, HybridSystems using LazySets using Reachability: solve using Plots using TaylorIntegration function bball_up!(t, x, dx) dx[1] = x[2] dx[2] = -9.8 - 0.1*(x[1])^2 return dx end 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") 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)

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

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)

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)

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)

## 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")

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)]

The first step is to normalize the polynomial (

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

Now collect the center and generators of the Zonotope (

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)

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)

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

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")

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)]

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]

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

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

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)

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)

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

### 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