Merge pull request #688 from JuliaReach/schillic/686

#686 - Decomposition for ASB07

src/ReachSets/ContinuousPost/ASB07_decomposed/ASB07_decomposed.jl

@@ -0,0 +1,22 @@
+export ASB07_decomposed
+
+struct ASB07_decomposed <: ContinuousPost
+    options::TwoLayerOptions
+
+    function ASB07_decomposed(𝑂::Options)
+        𝑂new = validate_and_wrap_options(𝑂, options_ASB07_decomposed())
+        return new(𝑂new)
+    end
+end
+
+# convenience constructor from pairs of symbols
+ASB07_decomposed(𝑂::Pair{Symbol,<:Any}...) =
+    ASB07_decomposed(Options(Dict{Symbol,Any}(𝑂)))
+
+# default options (they are added in the function validate_and_wrap_options)
+ASB07_decomposed() = ASB07_decomposed(Options())
+
+include("options.jl")
+include("init.jl")
+include("post.jl")
+include("reach.jl")

src/ReachSets/ContinuousPost/ASB07_decomposed/init.jl

@@ -0,0 +1,16 @@
+init(𝒜::ASB07_decomposed, 𝑆::AbstractSystem, 𝑂::Options) = init!(𝒜, 𝑆, copy(𝑂))
+
+function init!(𝒜::ASB07_decomposed, 𝑆::AbstractSystem, 𝑂::Options)
+    𝑂copy = copy(𝑂)
+    𝑂copy[:n] = statedim(𝑆)
+
+    𝑂validated = validate_solver_options_and_add_default_values!(𝑂copy)
+
+    𝑂validated[:partition] = 𝒜.options[:partition]
+
+    # compute blocks
+    𝑂validated[:blocks] =
+        compute_blocks(𝒜.options[:vars], 𝑂validated[:partition])
+
+    return 𝑂validated
+end

src/ReachSets/ContinuousPost/ASB07_decomposed/options.jl

@@ -0,0 +1,30 @@
+function options_ASB07_decomposed()
+    𝑂spec = Vector{OptionSpec}()
+
+    push!(𝑂spec, OptionSpec(:δ, 1e-2, domain=Float64, aliases=[:sampling_time],
+        domain_check=(v  ->  v > 0.), info="time step"))
+
+    push!(𝑂spec, OptionSpec(:max_order, 10, domain=Int,
+        info="maximum allowed order of zonotopes"))
+
+    push!(𝑂spec, OptionSpec(:order_discretization, 2, domain=Int,
+        info="order of Taylor approximation in discretization"))
+
+    push!(𝑂spec, OptionSpec(:set_operations_discretization, "zonotope",
+        domain=String, domain_check=(v  ->  v ∈ ["zonotope", "lazy"]),
+        info="type of set operations applied during discretization"))
+
+    push!(𝑂spec, OptionSpec(:block_options_init, Hyperrectangle, domain=Any,
+        info="option for the decomposition of the initial states"))
+
+    push!(𝑂spec, OptionSpec(:partition, [Int[]],
+        domain=AbstractVector{<:AbstractVector{Int}}, domain_check=ispartition,
+        info="block partition; a block is represented by a vector containing " *
+             "its indices"))
+
+    push!(𝑂spec, OptionSpec(:vars, Int[], domain=AbstractVector{Int},
+        domain_check=(v  ->  length(v) > 0 && all(e -> e > 0, v)),
+        info="variables of interest; default: all variables"))
+
+    return 𝑂spec
+end

src/ReachSets/ContinuousPost/ASB07_decomposed/post.jl

@@ -0,0 +1,56 @@
+function post(𝒜::ASB07_decomposed,
+              𝑃::InitialValueProblem{<:AbstractContinuousSystem},
+              𝑂::Options)
+    # =================================
+    # Initialization and discretization
+    # =================================
+
+    𝑂 = merge(𝒜.options.defaults, 𝑂, 𝒜.options.specified)
+    δ, T = 𝑂[:δ], 𝑂[:T]
+    N = round(Int, T / δ)
+    n = 𝑂[:n]
+    partition = 𝑂[:partition]
+    blocks = 𝑂[:blocks]
+    max_order = 𝑂[:max_order]
+
+    # compute and unrwap discretized system
+    𝑃_discrete = discretize(𝑃, δ; algorithm="interval_matrix",
+                            order=𝑂[:order_discretization],
+                            set_operations=𝑂[:set_operations_discretization])
+    Ω0, Φ = 𝑃_discrete.x0, 𝑃_discrete.s.A
+
+    # decompose initial states
+    # TODO this should use the concrete linear_map projection (see LazySets#1726)
+    Ωhat0 = array(decompose(Ω0, partition, 𝑂[:block_options_init]))
+
+    # ====================
+    # Flowpipe computation
+    # ====================
+
+    # preallocate output
+    Rsets = Vector{ReachSet{CartesianProductArray{Float64, LazySet{Float64}}}}(undef, N)
+
+    info("Reachable States Computation...")
+    @timing begin
+    if inputdim(𝑃_discrete) == 0
+        U = nothing
+    else
+        U = inputset(𝑃_discrete)
+    end
+    reach_ASB07_decomposed!(Rsets, Ωhat0, U, Φ, N, δ, max_order, n, partition,
+                            blocks)
+    end # timing
+
+    Rsol = ReachSolution(Rsets, 𝑂)
+
+    # ==========
+    # Projection
+    # ==========
+
+    if 𝑂[:project_reachset]
+        info("Projection...")
+        Rsol = @timing project(Rsol)
+    end
+
+    return Rsol
+end

src/ReachSets/ContinuousPost/ASB07_decomposed/reach.jl

@@ -0,0 +1,58 @@
+function reach_ASB07_decomposed!(R::Vector{<:ReachSet},
+                                 Ωhat0::Vector{<:LazySet},
+                                 U::Union{ConstantInput, Nothing},
+                                 ϕ::AbstractMatrix{IntervalArithmetic.Interval{NUM}},
+                                 N::Int,
+                                 δ::Float64,
+                                 max_order::Int,
+                                 n::Int,
+                                 partition,
+                                 blocks) where {NUM}
+    # initial reach set
+    t0, t1 = zero(δ), δ
+    R[1] = ReachSet(CartesianProductArray(Ωhat0[blocks]), t0, t1)
+
+    b = length(blocks)
+    Rₖ_array = Vector{LazySet{NUM}}(undef, b)
+    ϕpowerk = copy(ϕ)
+
+    if U != nothing
+        Whatk = Vector{Zonotope{NUM}}(undef, b)
+        inputs = next_set(U)
+        @inbounds for i in 1:b
+            bi = partition[blocks[i]]
+            Whatk[i] = linear_map(proj(bi, n), inputs)
+        end
+    end
+
+    k = 2
+    while k <= N
+        for i in 1:b
+            bi = partition[blocks[i]]
+            Rₖ_bi = ZeroSet(length(bi))
+            for (j, bj) in enumerate(partition)
+                block = IntervalMatrix(ϕpowerk[bi, bj])
+                if !iszero(block)
+                    Rₖ_bij = overapproximate(block * Ωhat0[j], Zonotope)
+                    Rₖ_bi = minkowski_sum(Rₖ_bi, Rₖ_bij)
+                end
+            end
+            if U != nothing
+                Rₖ_bi = minkowski_sum(Rₖ_bi, Whatk[i])
+            end
+            Rₖ_bi = reduce_order(Rₖ_bi, max_order)  # reduce order
+            Rₖ_array[i] = Rₖ_bi
+        end
+        Rₖ = CartesianProductArray(copy(Rₖ_array))
+
+        # store reach set
+        t0 = t1
+        t1 += δ
+        R[k] = ReachSet(Rₖ, t0, t1)
+
+        ϕpowerk *= ϕ
+
+        k += 1
+    end
+    return R
+end

src/ReachSets/ContinuousPost/BFFPSV18/reach_blocks.jl

@@ -29,14 +29,14 @@ given block indices.
 =#
 
 # helper functions
-@inline proj(bi::UnitRange{Int}, n::Int) =
+@inline proj(bi::AbstractVector{Int}, n::Int) =
         sparse(1:length(bi), bi, ones(length(bi)), length(bi), n)
 @inline proj(bi::Int, n::Int) = sparse([1], [bi], ones(1), 1, n)
-@inline row(ϕpowerk::AbstractMatrix, bi::UnitRange{Int}) = ϕpowerk[bi, :]
+@inline row(ϕpowerk::AbstractMatrix, bi::AbstractVector{Int}) = ϕpowerk[bi, :]
 @inline row(ϕpowerk::AbstractMatrix, bi::Int) = ϕpowerk[[bi], :]
-@inline row(ϕpowerk::SparseMatrixExp, bi::UnitRange{Int}) = get_rows(ϕpowerk, bi)
+@inline row(ϕpowerk::SparseMatrixExp, bi::AbstractVector{Int}) = get_rows(ϕpowerk, bi)
 @inline row(ϕpowerk::SparseMatrixExp, bi::Int) = Matrix(get_row(ϕpowerk, bi))
-@inline block(ϕpowerk_πbi::AbstractMatrix, bj::UnitRange{Int}) = ϕpowerk_πbi[:, bj]
+@inline block(ϕpowerk_πbi::AbstractMatrix, bj::AbstractVector{Int}) = ϕpowerk_πbi[:, bj]
 @inline block(ϕpowerk_πbi::AbstractMatrix, bj::Int) = ϕpowerk_πbi[:, [bj]]
 @inline store!(res, k, X, t0, t1, dimensions, N::Int) = (res[k] =
     SparseReachSet(X, t0, t1, dimensions))

src/ReachSets/ReachSets.jl

@@ -107,6 +107,8 @@ include("ContinuousPost/BFFPS19/BFFPS19.jl")
 
 include("ContinuousPost/ASB07/ASB07.jl")
 
+include("ContinuousPost/ASB07_decomposed/ASB07_decomposed.jl")
+
 # ========================
 # Reachability Algorithms
 # ========================

test/Reachability/solve_continuous.jl

@@ -221,9 +221,9 @@ property=(t,x)->x[2] <= 2.75
 𝑂 = Options(:T=>7.0, :mode=>"check", :property=>property)
 solve(𝑃, 𝑂, op=TMJets(:abs_tol=>1e-10, :orderT=>10, :orderQ=>2))
 
-# =====
-# ASB07
-# =====
+# ========================
+# ASB07 & ASB07_decomposed
+# ========================
 
 # example 1 from
 # "Reachability analysis of linear systems with uncertain parameters and inputs"
@@ -241,11 +241,21 @@ B = IntervalMatrix(hcat([1.0 ± 0.0; 1.0 ± 0.0]))
 U = ConstantInput(Interval(-0.05, 0.05))
 P_aff = IVP(CLCCS(A, B, nothing, U), X0)
 
+opC1 = ASB07()
+opC2 = ASB07(:δ => 0.04, :max_order => 10, :order_discretization => 4)
+opC3 = ASB07_decomposed(:vars => [1, 2], :partition => [[1], [2]])
+opC4 = ASB07_decomposed(:δ => 0.04, :max_order => 10,
+                        :order_discretization => 4, :vars => [1, 2],
+                        :partition => [[1], [2]])
+
 for P in [P_lin, P_aff]
     # default options
-    s = solve(P, Options(:T => 0.1), op=ASB07())
+    for opC in [opC1, opC3]
+        s = solve(P, Options(:T => 0.1), op=opC)
+    end
 
     # use specific options
-    opC = ASB07(:δ => 0.04, :max_order => 10, :order_discretization => 4)
-    s = solve(P, Options(:T => 5.0), op=opC)
+    for opC in [opC2, opC4]
+        s = solve(P, Options(:T => 5.0), op=opC)
+    end
 end