Best way to export results?

[willsharpless]

Hi all, thanks for all the help

With the previous L96 model (happy to share, omitted for brevity), I expanded the dimension to 40 and the analysis was able to complete at some point over night (wow, well done). Out of caution, I immediately exported the object returned by solve (using TMJets) and I used JLD, which I understand isn't the best but does the job.

It worked, however, the file is ~15 gb. Is there a better way to save the objects?

I understand this is use case dependent, but from the previous issues, maybe over-approximating with Zonotopes would yield a "lighter" structure and could be reformatted for future use if needed? While its nice to have the solver metadata, it's certainly an after thought.

[mforets]

It worked, however, the file is ~15 gb. Is there a better way to save the objects?

whoa, that's huge.

my first suggestion would indeed have been to export using JLD2.

I understand this is use case dependent, but from the previous issues, maybe over-approximating with Zonotopes would yield a "lighter" structure and could be reformatted for future use if needed?

if you are only interested in the plotting result and not in the taylor models themselves, then indeed i'd suggest first to convert to zonotope and then to save those.

in past projects we could handle lots (the order of a few millions) reach-sets, even for plotting. it's possible ,but i can be a bit tricky to do it efficiently.

[mforets]

(happy to share, omitted for brevity),

if it is only omitted for reason of brevity, then i encourage you to share it. we could take a look and potentially suggest optimizations.

[mforets]

I expanded the dimension to 40 and the analysis was able to complete at some point over night (wow, well done).

i thought the number of taylor series variables was bound to ~20, because of the size of the coeffs table

[willsharpless]

Here is the code! Annoyingly, the generated file can't be attached because jld is not a supported file type.

using LinearAlgebra, Interpolations, Plots, LaTeXStrings
using ReachabilityAnalysis, DifferentialEquations
gr()

const N = 40;
const nd = N;
const F = 8 * ones(N);

@taylorize function L96!(dx, x, p, t)

    ## L96 states
    dx[1] = (x[2] - x[N-1]) * x[N] - x[1] + F[1] + x[N+1]
    dx[2] = (x[3] - x[N])   * x[1] - x[2] + F[2] + x[N+2]
    for i = 3:(N-1)
        dx[i] = (x[i+1] - x[i-2]) * x[i-1] - x[i] + F[i] + x[N+i]
    end
    dx[N] = (x[1] - x[N-2]) * x[N-1] - x[N] + F[N] + x[N+N]

    ## Disturbance Dummy States
    local zerox = zero(x[1])
    for i in (N+1):(2N)
        dx[i] = zerox
    end

    return dx
end

### Observe Auto, No-D trajectory

t = 5.
xrand = vcat(rand(N), zeros(nd))
Xeq = solve(ODEProblem(L96!, xrand, (0.0, t)))[end][1:N]

## Differential Equations Solve
x0 = vcat(Xeq, zeros(nd))
prob = ODEProblem(L96!, x0, (0., t), saveat=0.:0.01:t)
sol = solve(prob); sola = hcat(sol.u...)
X = sola[1:N,:]

## Interpolate
itp = LinearInterpolation((1:N, sol.t), X)
xi_fine = range(extrema(1:N)..., length=300)
X_itp = [itp(xi,t) for t in sol.t, xi in xi_fine]

## Plot 
tfig = heatmap(xi_fine, sol.t, X_itp, cmap=cgrad(:RdBu), clims=(-11,11), ylims=[0, t], interpolate=true, title=latexstring("N=$N, F=$(F[1])"), dpi=300)
plot!(xlabel=L"x_i", ylabel=L"t\:\:(days)", xlims=[1, N], xticks=collect(1:2:N), yticks=(0:0.2:t, 0:Int(round(t/0.2))), legend=:topleft)

###  Taylor Model Solution with disturbance

t = 1.0

## Disturbance 
r𝒟 = 0.1;
c𝒟 = zeros(N);
Q𝒟 = inv(r𝒟) * diagm(ones(N));

X0 = Singleton(Xeq)
D0 = Hyperrectangle(; low=zeros(nd), high=zeros(nd))
D = Hyperrectangle(; low=c𝒟 - diag(inv(Q𝒟)), high=c𝒟 + diag(inv(Q𝒟)))

## Define Forwards Problem and Solve 

X̂0auto = X0 × D0;
X̂0 = X0 × D;
sys_auto = InitialValueProblem(BlackBoxContinuousSystem(L96!, length(low(X̂0auto))), X̂0auto);
sys = InitialValueProblem(BlackBoxContinuousSystem(L96!, length(low(X̂0))), X̂0);

@time sol_auto = solve(sys_auto; tspan=(0.0, t), alg=TMJets(; abstol=1e-12));
@time sol = solve(sys; tspan=(0.0, t), alg=TMJets(; abstol=1e-12));

solz_auto = overapproximate(sol_auto, Zonotope)
solz = overapproximate(sol, Zonotope)

## Plot

plot_vars = (1,2)
auto_plot = plot(sol_auto[end:-1:1]; vars=plot_vars, title=latexstring("Auto, t∈[0, $t]"))
dist_plot = plot(sol[end:-1:1]; vars=plot_vars, title=latexstring("|d| ≤ $r𝒟, t∈[0, $t]"));
v1 = plot(auto_plot, dist_plot, xlims=xlims(auto_plot), ylims=ylims(auto_plot));

plot_vars = (3,4)
auto_plot = plot(sol_auto[end:-1:1]; vars=plot_vars, title=latexstring("Auto, t∈[0, $t]"));
dist_plot = plot(sol[end:-1:1]; vars=plot_vars, title=latexstring("|d| ≤ $r𝒟, t∈[0, $t]"));
v2 = plot(auto_plot, dist_plot, xlims=xlims(auto_plot), ylims=ylims(auto_plot));

plot_vars = (5,6)
auto_plot = plot(sol_auto[end:-1:1]; vars=plot_vars, title=latexstring("Auto, t∈[0, $t]"));
dist_plot = plot(sol[end:-1:1]; vars=plot_vars, title=latexstring("|d| ≤ $r𝒟, t∈[0, $t]"));
v3 = plot(auto_plot, dist_plot, xlims=xlims(auto_plot), ylims=ylims(auto_plot));

ffig = plot(v1, v2, v3, layout=(3,1), size=(600, 800), dpi=300)

[willsharpless]

The ffig plot looks like this for the above,

[willsharpless]

I tried saving after over-approximating with Zonotope′s and it is better ... but still 14.28 gb lol

[schillic]

Wow, 80 dimensions. I am impressed that the method can handle that.

I tried saving after over-approximating with Zonotope′s and it is better ... but still 14.28 gb lol

How many sets are there? I expect for a time horizon t = 1.0 there should not be too many, and it looks like maybe 200 from the plots. So the size must come from the complexity of the individual sets. But still, several megabytes is unheard of...
Can you find out the order of one of the zonotopes (using order)?
You could try to overapproximate with a Hyperrectangle instead. That should really be cheap. But the plots will be less smooth.

[willsharpless]

Along the same lines, what is the most efficient way to access the vertices of any of these zonotopes? I'm doing this naively to compute the maximum norm in the set (I could also solve the convex opt problem, lmk if there is a better way built in).

Interestingly, the projection for the plots isn't too slow and ρ is fast for a given hyperplane, but vertices or vertices_list doesn't terminate.

Even when N=9 (making dim 18), vertices_list(set(solz[1])) doesn't terminate despite being a box around a Singleton, however, vertices_list(box_approximate(set(solz[1]))) then its very fast. The problem is the latter is often a big over approximation.

Part of the problem is this model form has all these dummy states (for disturbance) doubling the dimension, which do not change. Thus, I could project onto the first N dimensions but I'm not sure how to do that, nor if it will ultimately be enough.

[schillic]

@willsharpless: If you want the infinity norm, then the box approximation is what you want.

There is a function norm of a set, but it is rather naive and uses the vertices except for the infinity norm. So you can use that in your case.

[willsharpless]

Sorry missed this, meant 2-norm but eventually gave up on it and reduced to inf-norm