Using Taylor-model algorithm in a loop

[willsharpless]

Describe the bug

I would like to recursively solve and compute the hull of reach sets instead of one-shot solving with solve. However, right now when using TMJets (or not selecting the alg), the solution is outputted as a Flowpipe of TaylorModelReachSets and directly calling convexify() on this Flowpipe or an array of its elements, i.e. an array of TaylorModelReachSets, produces a no-method-matching error due the type. Is there a way to convert the TaylorModelReachSets in order to convexify() without overapproximation()?

I understand it is possible to overapproximate(sol, Zonotope) the solution, and then convexify, however as can be seen in the following example this often leads to large over approximations as the zonotope is (much) larger than the reach set.

To Reproduce

ReachabilityAnalysis v0.22.1

MWE: Note that the function VdP_ZT() recursively solves the backwards reachable set and computes the hull for some Δt step.

using LinearAlgebra, Plots
using ReachabilityAnalysis
plotlyjs()

function vanderpol!(dx, x, p, t)
    local μ = 1.0
    dx[1] = x[2]
    dx[2] = μ * (1.0 - x[1]^2) * x[2] + x[3]
    dx[3] = zero(x[3]) #input
    return dx
end

X0 = Hyperrectangle(; low=[1.25, 2.25], high=[1.55, 2.35])
U = Hyperrectangle(zeros(1), zeros(1)) 
X̂0 = X0 × U;

sys = @system(x' = vanderpol!(x), dim=3, x(0) ∈ X̂0)
sol = solve(sys, tspan=(1.0, 0.); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=6, orderQ=2)); #NSTEPS
# sol = solve(sys, tspan=(1.0, 0.)); # also produces TaylorModelReachSet
solz = overapproximate(sol, Zonotope);

# fig1 = plot(solz; vars=(1, 2), lw=0.3, title="Van der Pol", xlab="x₁", ylab="x₂");
# plot!(X0; vars=(1, 2), label="X0")

function VdP_ZT(X0, tf, Δt)

    Zi, ZTi = copy(X0), copy(X0);
    Zarr, ZTarr = [], [];

    for Δti = 1:Int(tf/Δt)

        ## Recursive Set Propagation
        sys = @system(x' = vanderpol!(x), dim=3, x(0) ∈ set(Zi));
        sol = solve(sys, tspan=(Δt, 0.); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=6, orderQ=2));
        solz = overapproximate(sol, Zonotope);
        Zi = solz[end];
        # Zi = sol[end];
        push!(Zarr, Zi); 

        ## Recursive Tube Propagation
        sys = @system(x' = vanderpol!(x), dim=3, x(0) ∈ set(ZTi));
        sol = solve(sys, tspan=(Δt, 0.); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=6, orderQ=2));
        solz = overapproximate(sol, Zonotope);
        ZTi = convexify([solz[1], solz[end]]);
        push!(ZTarr, ZTi)
    end

    return ZTi, Zarr, ZTarr
end

X̂0 = Hyperrectangle(; low=[1.25, 2.25, 0.0], high=[1.55, 2.35, 0.0])
tf, Δt = 1.0, 0.1
ZTi, Zarr, ZTarr = VdP_ZT(X̂0, tf, Δt); 

plot(ZTi; vars=(1,2), label="BRZT(T, tf)", lw=3.0, title="BRZ vs. BRZT: Δt=$Δt, tf=$tf", xlab="x¹₁", ylab="x¹₂");

for i=length(ZTarr)-1:-1:1
    label = i == length(ZTarr)-1 ? "BRZT(T, ti)" : ""
    plot!(ZTarr[i], vars=(1,2), label=label, color="orange", lw=3.0);
end

for i=length(ZTarr)-1:-1:1
    label = i == length(ZTarr)-1 ? "BRZ(T, ti)" : ""
    plot!(Zarr[i], vars=(1,2), label=label, color="green", lw=3.0);
end

plot!(solz, vars=(1,2), label="BRZ(T, ti)", color="green", lw=0.3)
plot!(X0, label="X0", color="magenta", lw=3)

Screenshots

When we compute the set recursively, we can compare what happens when we use overapproximate or not (top vs. bottom, in VdP_ZT(): Zi = solz[end]; vs. Zi = sol[end];). We observe that using overapproximate bloats the BRS compared to one-shot solving (which we don't want to do). In both cases, to convexify we have to over approximate and we can see this leads to a massive over approximation of the flowpipe.

Thank you for considering this! The package is a fantastic resource and I appreciate the hard work put into this.

[schillic]

Thanks for the detailed explanation. That made it very easy to understand.

First let me note that Taylor models are very precise set representations that are however hard to work with. That is why we often overapproximate with a zonotope or hyperrectangle and continue with that simpler set representation. But, as you noted, this inherently loses precision, and hence you do not want to continue from this overapproximation in a solve loop.

So the way to go is this: preserve the TaylorModelReachSet for the solve argument, and use the overapproximation for other calculations (like plotting, which implicitly also converts to zonotopes -- so what you see in the plots are not the actual Taylor models but overapproximations too).

The single call to solve at the beginning is doing the following under the hood: Since a TaylorModelReachSet is a set representation in time, one needs to evaluate that set at the final time point in order to get a precise (Taylor-model) set to continue with. But then you need to add time again to get the right type, so you basically have to unwrap and wrap again.

Coding this up is quite an adventure, so let me show you how to imitate¹ this (I modified your script):

using LinearAlgebra, Plots
using ReachabilityAnalysis

function vanderpol!(dx, x, p, t)
    local μ = 1.0
    dx[1] = x[2]
    dx[2] = μ * (1.0 - x[1]^2) * x[2] + x[3]
    dx[3] = zero(x[3]) #input
    return dx
end

X0 = Hyperrectangle(; low=[1.25, 2.25], high=[1.55, 2.35])
U = Hyperrectangle(zeros(1), zeros(1))
X̂0 = X0 × U;

sys = @system(x' = vanderpol!(x), dim=3, x(0) ∈ X̂0)
sol = solve(sys, tspan=(0.0, 1.); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=6, orderQ=2)); #NSTEPS
solz = overapproximate(sol, Zonotope);

fig1 = plot(solz; vars=(1, 2), lw=0.3, title="Van der Pol", xlab="x₁", ylab="x₂");
plot!(X0; vars=(1, 2), label="X0")

function VdP_ZT(X0, tf, Δt)

    Zi, ZTi = copy(X0), copy(X0);
    Zarr, ZTarr = [], [];
    # new
    TMi = copy(X0);
    ZTMarr = [];
    orderT = 6;
    n = dim(X0);
    zeroI = ReachabilityAnalysis.zeroI

    for Δti = 1:Int(tf/Δt)

        ## Recursive Set Propagation
        sys = @system(x' = vanderpol!(x), dim=3, x(0) ∈ set(Zi));
        sol = solve(sys, tspan=(0., Δt); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=6, orderQ=2));
        solz = overapproximate(sol, Zonotope);
        Zi = solz[end];
        # Zi = sol[end];
        push!(Zarr, Zi);

        ## Recursive Tube Propagation
        sys = @system(x' = vanderpol!(x), dim=3, x(0) ∈ set(ZTi));
        sol = solve(sys, tspan=(0., Δt); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=6, orderQ=2));
        solz = overapproximate(sol, Zonotope);
        ZTi = convexify([solz[1], solz[end]]);
        push!(ZTarr, ZTi)

        # new
        sys = @system(x' = vanderpol!(x), dim=3, x(0) ∈ TMi);
        sol = solve(sys, tspan=(0., Δt); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=orderT, orderQ=2));
        # last reach set (in time)
        last_tset = sol[end]
        # last reach set at last time point (without time)
        last_set = evaluate(set(last_tset), sup(domain(last_tset)))
        # convert polynomials back to a TaylorModelN vector (still without time)
        rem = remainder(last_tset)
        zB = ReachabilityAnalysis.zeroBox(n)
        sB = ReachabilityAnalysis.symBox(n)
        Y = [ReachabilityAnalysis.fp_rpa(ReachabilityAnalysis.TaylorModelN(last_set[j], rem[j], zB, sB)) for j in 1:n]
        # wrap TaylorModelN vector back into TaylorModel1 vector (in time)
        TMi = [ReachabilityAnalysis.TaylorModel1(Taylor1(polynomial(Y[j]), orderT), zeroI, zeroI, zeroI) for j in 1:n]
        solz = overapproximate(sol, Zonotope);
        ZTMi = convexify([solz[1], solz[end]]);
        push!(ZTMarr, ZTMi)
    end

    return ZTi, Zarr, ZTarr, ZTMarr
end

X̂0 = Hyperrectangle(; low=[1.25, 2.25, 0.0], high=[1.55, 2.35, 0.0])
tf, Δt = 1.0, 0.1
ZTi, Zarr, ZTarr, ZTMarr = VdP_ZT(X̂0, tf, Δt);

plot(ZTi; vars=(1,2), label="BRZT(T, tf)", lw=3.0, title="BRZ vs. BRZT: Δt=$Δt, tf=$tf", xlab="x¹₁", ylab="x¹₂")

for i=length(ZTarr)-1:-1:1
    label = i == length(ZTarr)-1 ? "BRZT(T, ti)" : ""
    plot!(ZTarr[i], vars=(1,2), label=label, color="orange", lw=3.0);
end

for i=length(ZTarr)-1:-1:1
    label = i == length(ZTarr)-1 ? "BRZ(T, ti)" : ""
    plot!(Zarr[i], vars=(1,2), label=label, color="green", lw=3.0);
end

for i=length(ZTMarr):-1:1
    label = i == length(ZTarr) ? "BRZTM(T, ti)" : ""
    plot!(ZTMarr[i], vars=(1,2), label=label, color="purple", lw=3.0);
end

# plot!(solz, vars=(1,2), label="BRZ(T, ti)", color="white", lw=0.3)
plot!(X0, label="X0", color="magenta", lw=3)

The plot shapes look different because your tspan definition was wrong: tspan=(Δt, 0.) should be tspan=(0., Δt) (lower bound, then upper bound). I opened an issue to catch this in the future.

¹There is one minor conceptual difference between a single call to solve and calls in a loop, but I do not expect any significant difference in practice.

[willsharpless]

Thank you a ton for this! I get what you mean. The hull you added to the plot looks much better, which I understand is basically the hull of the flowpipe now. Wrt to the tspan, I was actually solving the backwards zonotopes, hence the "BRZ"′s, and I was following the example in the docs, which doesn't seem to require a negative t for tspan=(t,0), but this doesn't change the solution you gave me so thank you nonetheless.

While this is in front of us, I admit that I was trying to do this to extra-safely, over-approximate the flowpipe, with the assumption that as Δt → 0, the hull of the Δt reach set and the initial set will contain all trajectories. I understand this is perhaps gratuitous as it appears the TaylorModelReachSet structure/alg accounts for the remainder and perhaps adds it as an extra disturbance to guarantee over-approximation, such as proposed in Althoff et al. 2008 or Koller et al. 2019 for example. However, the cites here do this with a linear model and yet the taylor series model used in the pkg is a high-order polynomial, no?

Do you think you could briefly explain this and/or provide some citations for the rigorous over-approx guarantees? The Taylor Method page looks like it will describe this explicitly but is (as any great, new codebase's manual is) still in development.

[schillic]

I was actually solving the backwards zonotopes, hence the "BRZ"′s, and I was following the example in the docs, which doesn't seem to require a negative t for tspan=(t,0)

I totally forgot this was allowed. When I played with your script, I believe something did not work because of that. So it may be that the TMJets algorithm does not support this feature properly (but I may be wrong). In any case, you can get a backward result for x' = f(x) simply by flipping the sign to x' = -f(x).

as Δt → 0, the hull of the Δt reach set and the initial set will contain all trajectories

This is actually guaranteed for any Δt.

However, the cites here do this with a linear model and yet the taylor series model used in the pkg is a high-order polynomial, no?

Exactly. The TMJets algorithm is especially good for nonlinear systems.

Do you think you could briefly explain this and/or provide some citations for the rigorous over-approx guarantees? The Taylor Method page looks like it will describe this explicitly but is (as any great, new codebase's manual is) still in development.

The Taylor method itself is a pretty standard approach to nonlinear reachability. We use the implementation in TaylorModels.jl. There you find the standard references in the README.

[lbenet]

The Taylor method itself is a pretty standard approach to nonlinear reachability. We use the implementation in TaylorModels.jl. There you find the standard references in the README.

Let me jump in to simply take full responsibility of my lack of time to carry on with the development. In any case, if you want to help, you are very welcome and I can guide you and provide all details.

Very very short: we aim at a validated integration using a high order polynomial in de independent variable (t), centered on an initial condition, and be useful on a domain D of the dependent variables. We parametrize the dependence on the dependent variables also with a polynomial on those variables (Taylor expansion), around the initial condition. Validated integration here means that we check an inclusion property that guarantees the unicity of the solution, including the bounds (Schauder’s fixed point theorem). All this is obtained by two Taylor integrations (of different order), the first to obtain the polynomial approximation, and the second to have an educated guess of the residual.

[willsharpless]

Validated integration here means that we check an inclusion property that guarantees the unicity of the solution, including the bounds (Schauder’s fixed point theorem). All this is obtained by two Taylor integrations (of different order), the first to obtain the polynomial approximation, and the second to have an educated guess of the residual.

Ah I think I see. So this uses the taylor approximated evolutions with a remainder to give differential inclusions? E.g. with some guarantees like in Scott & Barton 2013. Also is this estimation of the remainder a guaranteed upper-bound?

Sorry for the detailed questions, I was confused at first, thinking that this was using zonotopes to do the computations and then it must involve some lineariztion but I see now. Thanks for answering them, I appreciate it. I will use the standard citation for the package but please let me know if there is any papers you would like me to cite!

[mforets]

Sorry for the detailed questions, I was confused at first, thinking that this was using zonotopes to do the computations and then it must involve some lineariztion but I see now.

ReachabilityAnalysis.jl has a toy implementation using linearization in https://github.com/JuliaReach/ReachabilityAnalysis.jl/blob/master/src/Algorithms/QINT/reach_homog.jl but the n-dimensional one is not implemented. Though if worked from scratch, I'd recommend using a more SOTA approach like https://mediatum.ub.tum.de/doc/1620128/1620128.pdf

[mforets]

Why do you need to solve in a loop in the first place? That's not clear to me when reading the code.

If you do need to have different calls to solve with different initial sets, the key is in Christian's comments and code examples

So the way to go is this: preserve the TaylorModelReachSet for the solve argument, and use the overapproximation for other calculations (like plotting, which implicitly also converts to zonotopes -- so what you see in the plots are not the actual Taylor models but overapproximations too).

It is quite technical, but the best way to not loose in overapproximation is to try re-using the Taylor model in subsequent calls, instead of converting to zonotope and using that as initial set. We used such technique for https://ojs.aaai.org/index.php/AAAI/article/view/20790

[mforets]

Ah I think I see. So this uses the taylor approximated evolutions with a remainder to give differential inclusions? E.g. with some guarantees like in Scott & Barton 2013. Also is this estimation of the remainder a guaranteed upper-bound?

Yes, at least in theory TMJets already returns guaranteed estimates. In practice it is quite hard to implement, and we are aware of some bugs (I couldn't find the issue, not sure it's still open).

You can always compare with results obtained from another solver such as Flow* -- it is available in this library via FLOWSTAR, see https://github.com/JuliaReach/ReachabilityAnalysis.jl/tree/master/src/Algorithms/FLOWSTAR ).

Both are meant to perform rigorous integration, both in terms of remainder estimates and in terms of using interval arithmetic computations to handle floating point errors.