Aadesh Deshmukh / Jun 19 2019
Remix of

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

1. : 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,
2. 0 implies that when velocity becomes zero the ball starts to fall down.
```using Pkg
Pkg.clone("https://github.com/JuliaReach/Reachability.jl.git")

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:

1. the position is [4.9, 5.1]
2. 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
```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)```

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
using LazySets, Plots

# box (ball in the infinity norm)
B = BallInf([2.5,1.0], 0.5)
plot(B)
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

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

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

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