include compressed Kronecker form

src/Algorithms/CARLIN/CARLIN.jl

@@ -10,7 +10,7 @@ Implementation of the reachability method using Carleman linearization from [1].
 """
 Base.@kwdef struct CARLIN <: AbstractContinuousPost
     N::Int=2 # order of the algorithm
-    compressed::Bool=false # choose to use compressed Kronecker form
+    compress::Bool=false # choose to use compressed Kronecker form
     δ::Float64=0.1 # step size for the linear reachability solver
     bloat::Bool=true # choose to include the error estimate in the result
     resets::Union{Int,Vector{Float64}}=0 # choose the number of resets (equal spacing) or a vector specifying the reset times within the time interval [0, T]

src/Algorithms/CARLIN/kronecker.jl

@@ -4,6 +4,11 @@
 
 using LinearAlgebra: checksquare
 
+# TODO refactor to CarlemanLinearization.jl
+import CarlemanLinearization: lift_vector
+
+lift_vector(X0::IA.Interval, N) = lift_vector(Interval(X0), N)
+
 """
     kron_pow(x::IA.Interval, pow::Int)
 

src/Algorithms/CARLIN/post.jl

@@ -13,7 +13,12 @@ function _template(; n, N)
     return CustomDirections(dirs)
 end
 
-# TODO refactor
+# TODO check if it can be removed
+function _project(sol, vars)
+    πsol_1n = Flowpipe([ReachSet(set(project(R, vars)), tspan(R)) for R in sol])
+end
+
+# TODO refactor to MathematicalSystems
 import MathematicalSystems: statedim
 
 struct CanonicalQuadraticForm{T, MT<:AbstractMatrix{T}} <: AbstractContinuousSystem
@@ -30,12 +35,10 @@ function canonical_form(s::BlackBoxContinuousSystem)
     @requires Symbolics
     f = s.f
     # differentiate
-
 end
 
 function post(alg::CARLIN, ivp::IVP{<:AbstractContinuousSystem}, tspan; Δt0=interval(0), kwargs...)
-    @unpack N, compressed, δ, bloat, resets = alg
-    @assert compressed == false
+    @unpack N, compress, δ, bloat, resets = alg
 
     # for now we assume there are no resets
     if haskey(kwargs, :NSTEPS)
@@ -55,5 +58,5 @@ function post(alg::CARLIN, ivp::IVP{<:AbstractContinuousSystem}, tspan; Δt0=int
     dirs = _template(n=n, N=N)
     alg = LGG09(δ=δ, template=dirs, approx_model=Forward())
 
-    return reach_CARLIN_alg(X0, F1, F2; alg, resets, N, T, Δt0, bloat)
+    return reach_CARLIN_alg(X0, F1, F2; alg, resets, N, T, Δt0, bloat, compress)
 end

src/Algorithms/CARLIN/reach.jl

@@ -1,31 +1,31 @@
-function _project(sol, vars)
-    πsol_1n = Flowpipe([ReachSet(set(project(R, vars)), tspan(R)) for R in sol])
-end
-
 # general method given a reachability algorithm for the linear system
-function reach_CARLIN_alg(X0, F1, F2; alg, resets, N, T, Δt0, bloat)
+function reach_CARLIN_alg(X0, F1, F2; alg, resets, N, T, Δt0, bloat, compress)
     if resets == 0
-        reach_CARLIN(X0, F1, F2; alg, N, T, Δt=Δt0, bloat)
+        reach_CARLIN(X0, F1, F2; alg, N, T, Δt=Δt0, bloat, compress)
     else
-        reach_CARLIN_resets(X0, F1, F2, resets; alg, N, T, Δt=Δt0, bloat)
+        reach_CARLIN_resets(X0, F1, F2, resets; alg, N, T, Δt=Δt0, bloat, compress)
     end
 end
 
-function reach_CARLIN(X0, F1, F2; alg, N, T, Δt, bloat, A=nothing)
+function reach_CARLIN(X0, F1, F2; alg, N, T, Δt, bloat, compress, A=nothing)
     error_bound_func = error_bound_specabs
 
     # lift initial states
-    n = dim(X0)
-    ŷ0 = kron_pow_stack(X0, N) |> box_approximation
+    if compress
+        ŷ0 = lift_vector(X0, N) # see CarlemanLinearization/linearization.jl
+    else
+        ŷ0 = kron_pow_stack(X0, N) |> box_approximation
+    end
 
     # solve continuous ODE
     if isnothing(A)
-        A = build_matrix(F1, F2, N)
+        A = build_matrix(F1, F2, N; compress)
     end
     prob = @ivp(ŷ' = Aŷ, ŷ(0) ∈ ŷ0)
     sol = solve(prob, T=T, alg=alg, Δt0=Δt)
 
     # projection onto the first n variables
+    n = dim(X0)
     πsol_1n = _project(sol, 1:n)
 
     if !bloat
@@ -60,17 +60,17 @@ function _compute_resets(resets::Vector{Float64}, T)
     return [interval(aux[i], aux[i+1]) for i in 1:length(aux)-1]
 end
 
-function reach_CARLIN_resets(X0, F1, F2, resets; alg, N, T, Δt, bloat)
+function reach_CARLIN_resets(X0, F1, F2, resets; alg, N, T, Δt, bloat, compress)
 
     # build state matrix (remains unchanged upon resets)
-    A = build_matrix(F1, F2, N)
+    A = build_matrix(F1, F2, N; compress)
 
     # time intervals to compute
     time_intervals = _compute_resets(resets, T)
 
     # compute until first chunk
     T1 = sup(first(time_intervals))
-    sol_1 = reach_CARLIN(X0, F1, F2; alg, N, T=T1, Δt=interval(0)+Δt, bloat, A=A)
+    sol_1 = reach_CARLIN(X0, F1, F2; alg, N, T=T1, Δt=interval(0)+Δt, bloat, compress, A=A)
 
     # preallocate output flowpipe
     fp_1 = flowpipe(sol_1)
@@ -85,7 +85,7 @@ function reach_CARLIN_resets(X0, F1, F2, resets; alg, N, T, Δt, bloat)
     for i in 2:length(time_intervals)
         T0 = T1
         Ti = sup(time_intervals[i])
-        sol_i = reach_CARLIN(X0, F1, F2; alg, N, T=Ti-T0, Δt=interval(T0)+Δt, A=A, bloat)
+        sol_i = reach_CARLIN(X0, F1, F2; alg, N, T=Ti-T0, Δt=interval(T0)+Δt, bloat, compress, A=A)
         push!(out, flowpipe(sol_i))
         X0 = box_approximation(set(sol_i[end]))
         T1 = Ti

test/algorithms/CARLIN.jl

@@ -13,4 +13,9 @@
     sol = solve(prob, T=10.0, alg=CARLIN(N=3, bloat=false))
     @test dim(sol) == 1
     @test ρ([1.0], sol(10.0)) < 0.01
+
+    # Use the compressed Kronecker form
+    sol = solve(prob, T=10.0, alg=CARLIN(N=3, bloat=false, compress=true))
+    @test dim(sol) == 1
+    @test ρ([1.0], sol(10.0)) < 0.01
 end