Merge pull request #626 from JuliaReach/carlin

Add initial version of the CARLIN algorithm

Project.toml

@@ -3,6 +3,7 @@ uuid = "1e97bd63-91d1-579d-8e8d-501d2b57c93f"
 version = "0.18.0"
 
 [deps]
+CarlemanLinearization = "4803f6b2-022a-4c1b-a771-522a3413ec86"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 ExprTools = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
 HybridSystems = "2207ec0c-686c-5054-b4d2-543502888820"
@@ -23,6 +24,7 @@ TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 
 [compat]
+CarlemanLinearization = "0.3"
 CDDLib = "0.5 - 0.9"
 CommonSolve = "0.2"
 DifferentialEquations = "6.17, 7"

src/Algorithms/CARLIN/CARLIN.jl

@@ -1,5 +1,21 @@
-# refactored to JuliaReach/CarlemanLinearization.jl
-# include("error_bounds.jl")
+using CarlemanLinearization
+
+"""
+    CARLIN <: AbstractContinuousPost
+
+Implementation of the reachability method using Carleman linearization from [1].
+
+[1] Forets, Marcelo, and Christian Schilling. "Reachability of weakly nonlinear systems using Carleman linearization."
+    International Conference on Reachability Problems. Springer, Cham, 2021.
+"""
+Base.@kwdef struct CARLIN <: AbstractContinuousPost
+    N::Int=2 # order of the algorithm
+    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]
+end
 
-# extends functionality in JuliaReach/CarlemanLinearization.jl to work with LazySets.jl types
 include("kronecker.jl")
+include("post.jl")
+include("reach.jl")

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

@@ -0,0 +1,62 @@
+function _template(; n, N)
+    dirs = Vector{SingleEntryVector{Float64}}()
+    d = sum(n^i for i in 1:N)
+
+    for i in 1:n
+        x = SingleEntryVector(i, d, 1.0)
+        push!(dirs, x)
+    end
+    for i in 1:n
+        x = SingleEntryVector(i, d, -1.0)
+        push!(dirs, x)
+    end
+    return CustomDirections(dirs)
+end
+
+# 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
+    F1::MT
+    F2::MT
+end
+statedim(f::CanonicalQuadraticForm) = size(f.F1, 1)
+canonical_form(s::CanonicalQuadraticForm) = s.F1, s.F2
+export CanonicalQuadraticForm
+#-
+
+# TODO generalize to AbstractContinuousSystem using vector_field
+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, compress, δ, bloat, resets = alg
+
+    # for now we assume there are no resets
+    if haskey(kwargs, :NSTEPS)
+        NSTEPS = kwargs[:NSTEPS]
+        T = NSTEPS * δ
+    else
+        # get time horizon from the time span imposing that it is of the form (0, T)
+        T = _get_T(tspan, check_zero=true, check_positive=true)
+        NSTEPS = ceil(Int, T / δ)
+    end
+
+    # extract initial states
+    X0 = initial_state(ivp)
+    # compute canonical quadratic form from the given system
+    F1, F2 = canonical_form(system(ivp))
+    n = statedim(ivp)
+    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, compress)
+end

src/Algorithms/CARLIN/reach.jl

@@ -0,0 +1,95 @@
+# general method given a reachability algorithm for the linear system
+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, compress)
+    else
+        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, compress, A=nothing)
+    error_bound_func = error_bound_specabs
+
+    # lift initial states
+    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; 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
+        return πsol_1n
+    end
+
+    # compute errors
+    errfunc = error_bound_func(X0, Matrix(F1), Matrix(F2), N=N)
+
+    # evaluate error bounds for each reach-set in the solution
+    E = [errfunc.(tspan(R)) for R in sol]
+
+    # if the interval is always > 0 then we can just take max(Ei)
+
+    # symmetrize intervals
+    E_rad = [symmetric_interval_hull(Interval(ei)) for ei in E]
+    E_ball = [BallInf(zeros(n), max(ei)) for ei in E_rad]
+
+    # sum the solution with the error
+    fp_bloated = [ReachSet(set(Ri) ⊕ Ei, tspan(Ri)) for (Ri, Ei) in zip(πsol_1n, E_ball)] |> Flowpipe
+
+    return fp_bloated
+end
+
+function _compute_resets(resets::Int, T)
+    return mince(0 .. T, resets+1)
+end
+
+function _compute_resets(resets::Vector{Float64}, T)
+    # assumes initial time is 0
+    aux = vcat(0.0, resets, 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, compress)
+
+    # build state matrix (remains unchanged upon resets)
+    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, compress, A=A)
+
+    # preallocate output flowpipe
+    fp_1 = flowpipe(sol_1)
+    out = Vector{typeof(fp_1)}()
+    push!(out, fp_1)
+
+    # approximate final reach-set
+    Rlast = sol_1[end]
+    X0 = box_approximation(set(Rlast))
+
+    # compute remaining chunks
+    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, bloat, compress, A=A)
+        push!(out, flowpipe(sol_i))
+        X0 = box_approximation(set(sol_i[end]))
+        T1 = Ti
+    end
+
+    return MixedFlowpipe(out)
+end

src/Initialization/exports.jl

@@ -11,6 +11,7 @@ export
     ASB07,
     BFFPSV18,
     BOX,
+    CARLIN,
     GLGM06,
     INT,
     LGG09,

test/algorithms/CARLIN.jl

@@ -0,0 +1,21 @@
+@testset "CARLIN algorithm: logistic model" begin
+    r = -0.5
+    K = 0.8
+    F1 = hcat(r)
+    F2 = hcat(-r/K)
+    X0 = 0.47 .. 0.53
+    prob = @ivp(CanonicalQuadraticForm(F1, F2), x(0) ∈ X0)
+
+    sol = solve(prob, T=10.0, alg=CARLIN(N=3, bloat=true))
+    @test dim(sol) == 1
+    @test ρ([1.0], sol(10.0)) < 0.2
+
+    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

test/runtests.jl

@@ -41,6 +41,7 @@ include("symbolics.jl")
 
 include("algorithms/INT.jl")
 include("algorithms/BOX.jl")
+include("algorithms/CARLIN.jl")
 include("algorithms/GLGM06.jl")
 include("algorithms/LGG09.jl")
 include("algorithms/ASB07.jl")