Discrete sequence using Taylor models

discrete_seq.jl

using ReachabilityAnalysis
using ReachabilityAnalysis.TaylorModels
using ReachabilityAnalysis.TaylorSeries
using Random
using StatsBase

# -----------------
# Numbers
# -----------------

function run(::Type{Float64})
    # Choose scenario
    n = 5
    u = rand(range(0.2, 0.5, length=10), n)
    x0 = rand(-1:0.01:1 * 0.01)

    # Propagate
    x = zeros(n+1)
    x[1] = x0
    for k in 1:n
        x[k+1] = x[k]^2 + x[k] + u[k]
    end
    x
end
run() = run(Float64)

out = [run()[5] for _ in 1:20]
target = 3 .. 5
@show countmap(out .∈ target)

# ----------------------------------
# Taylor models, no uncertainty
# ----------------------------------

function run(::Type{TaylorModel1})
    # Choose scenario
    TaylorSeries.set_taylor1_varname("x")
    n = 5
    ord = 6

    dom = interval(-1, 1)
    u = rand(range(0.2, 0.5, length=10), n)
    U = [TaylorModel1(ui + 0.0 * Taylor1(Float64, ord), zero(dom), zero(dom), dom) for ui in u]
    x0 = rand(-1:0.01:1 * 0.01)    
    X0 = TaylorModel1(x0 + 0.0 * Taylor1(Float64, ord), zero(dom), zero(dom), dom)

    # Propagate
    X = Vector{TaylorModel1{Float64, Float64}}(undef, n+1)
    X[1] = X0
    for k in 1:n
        X[k+1] = X[k]^2 + X[k] + U[k]
    end
    X
end

out = [run(TaylorModel1)[5] for _ in 1:20]
target = 3 .. 5
ovalues = [evaluate(o, -1 .. 1) for o in out]
@show countmap(ovalues .⊆ target)


# -------------------------------------
# Taylor models with uncertainty in u
# -------------------------------------

function run(::Type{TaylorModel1}, W)
    # Choose scenario
    TaylorSeries.set_taylor1_varname("x")
    n = 5
    ord = 6

    dom = interval(-1, 1)
    u = rand(range(0.2, 0.5, length=10), n) .+ W
    x = Taylor1(Float64, ord)
    U = [TaylorModel1(ui.hi * x + ui.lo * (1 - x), zero(dom), zero(dom), dom) for ui in u]
    x0 = rand(-1:0.01:1 * 0.01)    
    X0 = TaylorModel1(x0 + 0.0 * x, zero(dom), zero(dom), dom)

    # Propagate
    X = Vector{TaylorModel1{Float64, Float64}}(undef, n+1)
    X[1] = X0
    for k in 1:n
        X[k+1] = X[k]^2 + X[k] + U[k]
    end
    X
end

W = interval(-0.005, 0.005)
out = [run(TaylorModel1, W)[5] for _ in 1:20]
target = 3 .. 5
ovalues = [evaluate(o, -1 .. 1) for o in out]
@show countmap(ovalues .⊆ target)