Polytopic overappoximation of a taylor model

[mforets]

Let F = (p(x₀, t), I) be an n-dimensional TM flowpipe, x₀ ∈ X₀, t ∈ [0, δ]. Given a set of template directions d1, ..., dr in Rn, we can evaluate the support functions

bi = sup{di^T x} s.t. x = p(x₀, t) + y , x₀ ∈ X₀ , t ∈ [0, δ], y ∈ I
for each i = 1, ..., r.

Note that one can find a value which is larger than bi using techniques from SDP relaxation (*).

Refs:

• Section 4.3, Chen, X. (2015). Reachability analysis of non-linear hybrid systems using taylor models (Doctoral dissertation, Fachgruppe Informatik, RWTH Aachen University)..
• (*) For a fixed relaxation degree one obtains an overapproximation of the true solution, which is good in this case, since it will give an enclosing polytope.
• See this notebook for an example.

[dfcaporale]

Our try with @mforets today:

using TaylorModels
const IA = IntervalArithmetic;

I = IA.Interval(-0.5, 0.5) # interval remainder

x₀ = IA.Interval(0.0) # expansion point
D = IA.Interval(-3.0, 1.0)
p1 = Taylor1([2.0, 1.0], 2) # define a linear polynomial
p2 = Taylor1([0.9, 3.0], 2) # define another linear polynomial
vTM = [TaylorModel1(pi, I, x₀, D) for pi in [p1, p2]] # define vector of Taylor models

using LazySets.Approximations:AbstractDirections

function _overapproximate(vTM::Vector{TaylorModel1{T, S}}, dir::AbstractDirections{T}) where {T, S}

    Y2 = Vector{HalfSpace{T, Vector{T}}}();
    m = length(vTM)

    model = Model(Ipopt.Optimizer)
    register(model, :dot, m, dot; autodiff = true)

    dom_vTM = domain(vTM[1]) # asumimos mismo dominio de t en todas las variables
    t = @variable(model, dom_vTM.lo ≤ t ≤ dom_vTM.hi) # variables de estado
    r = @variable(model, r[1:m]) # remainders
    for i = 1:length(vTM)
        set_upper_bound(r[i], remainder(vTM[i]).hi)
        set_lower_bound(r[i], remainder(vTM[i]).lo)
    end    
    
    for d in dirs
        y = dot(d,vTM)
        coefs = y.pol.coeffs
        @NLobjective(model, Max, coefs[1] + sum(coefs[i]*t^(i-1) for i in 2:length(coefs)) + 
                                 sum(r[i]*d[i] for i in 1:m))
        optimize!(model)
        o = objective_value(model)
        push!(Y2, HalfSpace(d, o))
    end
    
    return HPolyhedron(Y2)
end