From 6d5e645e5b10b23b6427b42e317ac256049b7242 Mon Sep 17 00:00:00 2001 From: Kostiantyn Potomkin Date: Wed, 26 Jun 2019 21:24:42 +1000 Subject: [PATCH] Added BFFPS19 and DecomposedDiscrete Post --- .../ContinuousPost/BFFPS19/BFFPS19.jl | 291 +++++++++++++++++ src/ReachSets/ContinuousPost/BFFPS19/reach.jl | 296 +++++++++++++++++ .../ContinuousPost/BFFPS19/reach_blocks.jl | 301 ++++++++++++++++++ .../DiscretePost/DecomposedDiscretePost.jl | 144 +++++++++ src/ReachSets/ReachSets.jl | 3 + src/solve.jl | 12 +- 6 files changed, 1045 insertions(+), 2 deletions(-) create mode 100644 src/ReachSets/ContinuousPost/BFFPS19/BFFPS19.jl create mode 100644 src/ReachSets/ContinuousPost/BFFPS19/reach.jl create mode 100644 src/ReachSets/ContinuousPost/BFFPS19/reach_blocks.jl create mode 100644 src/ReachSets/DiscretePost/DecomposedDiscretePost.jl diff --git a/src/ReachSets/ContinuousPost/BFFPS19/BFFPS19.jl b/src/ReachSets/ContinuousPost/BFFPS19/BFFPS19.jl new file mode 100644 index 00000000..8cc4f92c --- /dev/null +++ b/src/ReachSets/ContinuousPost/BFFPS19/BFFPS19.jl @@ -0,0 +1,291 @@ +export BFFPS19 + +# =============================================================== +# Bogomolov, Forets, Frehse, Potomkin, Schilling. +# =============================================================== + +function options_BFFPS19() + return OptionSpec[ + # general options + OptionSpec(:discretization, "forward", domain=String, + info="model for bloating/continuous time analysis"), + OptionSpec(:δ, 1e-2, domain=Float64, aliases=[:sampling_time], + domain_check=(v -> v > 0.), info="time step"), + OptionSpec(:vars, Int[], domain=AbstractVector{Int}, domain_check=( + v -> length(v) > 0 && all(e -> e > 0, v)), + info="variables of interest; default: all variables"), + OptionSpec(:guards_proj, Vector(), domain=AbstractVector, + info="internal storage of projected guard constraints"), + OptionSpec(:partition, [Int[]], + domain=AbstractVector{<:AbstractVector{Int}}, domain_check= + ispartition, + info="block partition; a block is represented by a vector " * + "containing its indices"), + + # discretization options + OptionSpec(:sih_method, "concrete", domain=String, + info="method to compute the symmetric interval hull in discretization"), + OptionSpec(:exp_method, "base", domain=String, + info="method to compute the matrix exponential"), + OptionSpec(:assume_sparse, false, domain=Bool, + info="use an analysis for sparse discretized matrices?"), + + # reachability options + OptionSpec(:lazy_inputs_interval, lazy_inputs_interval_always, + domain=Union{Int, Function}, + domain_check=(v -> !(v isa Int) || v >= -1), + info="length of interval in which the inputs are handled as a " * + "lazy set (``-1`` for 'never'); may generally also be a " * + "predicate over indices; the default corresponds to ``-1``"), + + # approximation options + OptionSpec(:block_options, nothing, domain=Any, + info="short hand to set ':block_options_init' and " * + "':block_options_iter'"), + OptionSpec(:block_options_init, nothing, domain=Any, + info="option for the approximation of the initial states " * + "(during decomposition)"), + OptionSpec(:block_options_iter, nothing, domain=Any, + info="option for the approximation of the states ``X_k``, ``k>0``"), + + # convenience options + OptionSpec(:assume_homogeneous, false, domain=Bool, + info="ignore dynamic inputs during the analysis?"), + ] +end + +function normalization_BFFPS19!(𝑂::TwoLayerOptions) + # :lazy_inputs_interval option: convert integers to functions + if haskey_specified(𝑂, :lazy_inputs_interval) + v = 𝑂[:lazy_inputs_interval] + if v isa Int + if v == -1 + 𝑂.specified[:lazy_inputs_interval] = lazy_inputs_interval_never + elseif v == 0 + 𝑂.specified[:lazy_inputs_interval] = lazy_inputs_interval_always + else + 𝑂.specified[:lazy_inputs_interval] = (k -> k % v == 0) + end + end + end + + # :block_options options: use convenience option for '_init' and '_iter' + if !haskey_specified(𝑂, :block_options_init) && + haskey_specified(𝑂, :block_options) + 𝑂.specified[:block_options_init] = 𝑂[:block_options] + end + if !haskey_specified(𝑂, :block_options_iter) && + haskey_specified(𝑂, :block_options) + 𝑂.specified[:block_options_iter] = 𝑂[:block_options] + end + + nothing +end + +function validation_BFFPS19(𝑂) + # block_options + for b_options in [:block_options, :block_options_init, :block_options_iter] + if haskey_specified(𝑂, b_options) + bo = 𝑂[b_options] + if bo isa Type{<:LazySet} || bo isa Type{<:AbstractDirections} + # uniform options + elseif bo isa Symbol + # template directions + option = get(Utils.template_direction_symbols, bo, nothing) + if option == nothing + throw(DomainError(bo, "if the `$b_options` option " * + "is a Symbol, it must be one of " * + "$(keys(Utils.template_direction_symbols))")) + end + 𝑂.specified[b_options] = option + elseif bo isa AbstractVector || bo isa Dict{Int, Any} + # mapping + elseif bo isa Real || bo isa Pair{<:UnionAll, <:Real} + ε = bo isa Real ? bo : bo[2] + if ε <= 0 + throw(DomainError(ε, "the `$b_options` option must be " * + "positive")) + end + elseif b_options == :block_options_iter && bo == nothing + # no overapproximation + else + throw(DomainError(bo == nothing ? "nothing" : bo, + "the `$b_options` option does not accept the given input")) + end + end + end + + nothing +end + +""" + BFFPS19 <: ContinuousPost + +Implementation of the reachability algorithm for purely continuous linear +time-invariant systems using block decompositons by S. Bogomolov, M. Forets, +G. Frehse, A. Podelski, C. Schilling and F. Viry [1]. + +### Fields + +- `options` -- an `Options` structure that holds the algorithm-specific options + +### Notes + +The following options are available: + +```julia +$(print_option_spec(options_BFFPS19())) +``` + +### Algorithm + +We refer to [1] for technical details. + +[1] [Reach Set Approximation through Decomposition with Low-dimensional Sets +and High-dimensional Matrices](https://dl.acm.org/citation.cfm?id=3178128). +S. Bogomolov, M. Forets, G. Frehse, A. Podelski, C. Schilling, F. Viry. +HSCC '18 Proceedings of the 21st International Conference on Hybrid Systems: +Computation and Control (part of CPS Week). +""" +struct BFFPS19 <: ContinuousPost + options::TwoLayerOptions + + function BFFPS19(𝑂::Options) + normalized_𝑂 = validate_and_wrap_options(𝑂, options_BFFPS19(); + validation=validation_BFFPS19, + normalization=normalization_BFFPS19!) + return new(normalized_𝑂) + end +end + +# convenience constructor from pairs of symbols +BFFPS19(𝑂::Pair{Symbol,<:Any}...) = BFFPS19(Options(Dict{Symbol,Any}(𝑂))) + +# default options +BFFPS19() = BFFPS19(Options()) + +include("reach.jl") +include("reach_blocks.jl") + +init(𝒫::BFFPS19, 𝑆::AbstractSystem, 𝑂::Options) = init!(𝒫, 𝑆, copy(𝑂)) + +function init!(𝒫::BFFPS19, 𝑆::AbstractSystem, 𝑂::Options) + # state dimension for (purely continuous or purely discrete systems) + 𝑂copy = copy(𝑂) + 𝑂copy[:n] = statedim(𝑆) + + # solver-specific options (adds default values for unspecified options) + 𝑂validated = validate_solver_options_and_add_default_values!(𝑂copy) + + # :vars option; default: all variables + if haskey_specified(𝒫.options, :vars) + 𝑂validated[:vars] = 𝒫.options[:vars] + else + 𝑂validated[:vars] = 1:𝑂validated[:n] + end + + # :partition option: use 1D blocks + if haskey_specified(𝒫.options, :partition) + 𝑂validated[:partition] = 𝒫.options[:partition] + else + 𝑂validated[:partition] = [[i] for i in 1:𝑂validated[:n]] + end + + opD = 𝒫.options[:opD] + if opD isa DecomposedDiscretePost + HS = 𝒫.options[:HS] + constrained_dims = constrained_dimensions(HS) + out_vars = opD.options[:out_vars] + loc_id = 𝒫.options[:loc_id] + temp_vars = unique([out_vars; constrained_dims[loc_id]]) + opD.options[:temp_vars] = temp_vars + guards_constraints = [guard(HS, trans) for trans in out_transitions(HS, loc_id)] + 𝑂validated[:vars] = temp_vars + 𝑂validated[:guards_proj] = [project(c, temp_vars) for c in guards_constraints] + 𝑂validated[:blocks] = compute_blocks(𝑂validated[:vars], 𝑂validated[:partition]) + + info("- Total") + else + throw(MethodError("This continuous post operator works only with DecomposedDiscretePost")) + end + + # :block_options_init & :block_options_iter options: + # set default according to :partition + if !haskey_specified(𝒫.options, :block_options_init) + 𝑂validated[:block_options_init] = + compute_default_block_options(𝑂validated[:partition]) + end + if !haskey_specified(𝒫.options, :block_options_iter) + 𝑂validated[:block_options_iter] = + compute_default_block_options(𝑂validated[:partition]) + end + + # Input -> Output variable mapping + 𝑂validated[:inout_map] = inout_map_reach(𝑂validated[:partition], 𝑂validated[:blocks], 𝑂validated[:n]) + + if 𝑂validated[:project_reachset] + 𝑂validated[:output_function] = nothing + else + 𝑂validated[:output_function] = 𝑂validated[:projection_matrix] + end + + return 𝑂validated +end + +""" + post(𝒫::BFFPS19, 𝑆::AbstractSystem, 𝑂::Options) + +Calculate the reachable states of the given initial value problem using `BFFPS19`. + +### Input + +- `𝒫` -- post operator of type `BFFPS19` +- `𝑆` -- sytem, initial value problem for a continuous ODE +- `𝑂` -- algorithm-specific options +""" +function post(𝒫::BFFPS19, 𝑆::AbstractSystem, 𝑂_input::Options) + 𝑂 = TwoLayerOptions(merge(𝑂_input, 𝒫.options.specified), 𝒫.options.defaults) + + if 𝑂[:mode] == "reach" + info("Reachable States Computation...") + @timing begin + Rsets = reach_mixed(𝑆, 𝑂) + end + + # Projection + if 𝑂[:project_reachset] + info("Projection...") + RsetsProj = @timing project(Rsets, 𝑂) + else + RsetsProj = Rsets + end + + return ReachSolution(RsetsProj, 𝑂_input) + + elseif 𝑂[:mode] == "check" + info("invariants are currently not supported in 'check' mode") + + # Input -> Output variable mapping in property + property = inout_map_property(𝑂[:property], 𝑂[:partition], 𝑂[:blocks], 𝑂[:n]) + + # ================= + # Property checking + # ================= + info("Property Checking...") + @timing begin + answer = check_property(𝑆, property, 𝑂) + info("- Total") + end + + if answer == 0 + info("The property is satisfied!") + return CheckSolution(true, -1, 𝑂_input) + else + info("The property may be violated at index $answer," * + " (time point $(answer * 𝑂[:δ]))!") + return CheckSolution(false, answer, 𝑂_input) + end + else + error("unsupported mode $(𝑂[:mode])") + end # mode +end diff --git a/src/ReachSets/ContinuousPost/BFFPS19/reach.jl b/src/ReachSets/ContinuousPost/BFFPS19/reach.jl new file mode 100644 index 00000000..61f6fa80 --- /dev/null +++ b/src/ReachSets/ContinuousPost/BFFPS19/reach.jl @@ -0,0 +1,296 @@ +using ..Utils: LDS, CLCDS + +using LazySets: CacheMinkowskiSum, + isdisjoint + +import LazySets.Approximations: overapproximate +import ProgressMeter: update! + +""" +combine_cpas(cpa1::CartesianProductArray{N, S}, cpa2::CartesianProductArray{N, S}, + blocks1::Vector{Int}, blocks2::Vector{Int}) where {N, S<:LazySet{N}} + +Compose high-dimensional cartesian product array from 2 low-dimensional CPA's + +### Input + +- `cpa1` -- Cartesian Product Array +- `cpa2` -- Cartesian Product Array +- `blocks1` -- Block structure of cpa1 in sense of high-dimensional set +- `blocks2` -- Block structure of cpa2 in sense of high-dimensional set + +### Output + +A sequence of [`ReachSet`](@ref)s. + +### Notes + +A dictionary with available algorithms is available via +`Reachability.available_algorithms`. + +The numeric type of the system's coefficients and the set of initial states +is inferred from the first parameter of the system (resp. lazy set), ie. +`NUM = first(typeof(problem.s).parameters)`. +""" +function combine_cpas(cpa1::CartesianProductArray{N, S}, cpa2::CartesianProductArray{N, S}, + blocks1::Vector{Int}, blocks2::Vector{Int}) where {N, S<:LazySet{N}} + result = Vector{S}(undef, length(array(cpa1)) + length(array(cpa2))) + result[blocks1] = array(cpa1) + result[blocks2] = array(cpa2) + return CartesianProductArray(result) +end + +""" + reach(problem, options) + +Interface to reachability algorithms for an LTI system. + +### Input + +- `problem` -- initial value problem +- `options` -- additional options + +### Output + +A sequence of [`ReachSet`](@ref)s. + +### Notes + +A dictionary with available algorithms is available via +`Reachability.available_algorithms`. + +The numeric type of the system's coefficients and the set of initial states +is inferred from the first parameter of the system (resp. lazy set), ie. +`NUM = first(typeof(problem.s).parameters)`. +""" +function reach_mixed(problem::Union{IVP{<:CLDS{NUM}, <:LazySet{NUM}}, + IVP{<:CLCDS{NUM}, <:LazySet{NUM}}}, + options::TwoLayerOptions + )::Vector{<:ReachSet} where {NUM <: Real} + # optional matrix conversion + problem = matrix_conversion(problem, options) + + # list containing the arguments passed to any reachability function + args = [] + + # coefficients matrix + A = problem.s.A + push!(args, A) + + # determine analysis mode (sparse/dense) for lazy_expm mode + if A isa SparseMatrixExp + push!(args, Val(options[:assume_sparse])) + end + + n = statedim(problem) + blocks = options[:blocks] + partition = convert_partition(options[:partition]) + T = options[:T] + N = (T == Inf) ? nothing : ceil(Int, T / options[:δ]) + # Cartesian decomposition of the initial set + if length(partition) == 1 && length(partition[1]) == n && + options[:block_options_init] == LinearMap + info("- Skipping decomposition of X0") + Xhat0 = LazySet{NUM}[problem.x0] + else + info("- Decomposing X0") + @timing begin + Xhat0 = array(decompose(problem.x0, options[:partition], + options[:block_options_init])) + end + end + + # determine output function: linear map with the given matrix + output_function = options[:output_function] != nothing ? + (x -> options[:output_function] * x) : + nothing + + # preallocate output vector + if output_function == nothing + res_type = ReachSet{CartesianProductArray{NUM, LazySet{NUM}}, NUM} + else + res_type = ReachSet{Hyperrectangle{NUM}, NUM} + end + res = (N == nothing) ? Vector{res_type}() : Vector{res_type}(undef, N) + + # shortcut if only the initial set is required + if N == 1 + res[1] = res_type( + CartesianProductArray{NUM, LazySet{NUM}}(Xhat0[blocks]), + zero(NUM), options[:δ]) + return res + end + push!(args, Xhat0) + + # inputs + if !options[:assume_homogeneous] && inputdim(problem) > 0 + U = inputset(problem) + else + U = nothing + end + push!(args, U) + + # overapproximation function for states + block_options_iter = options[:block_options_iter] + if block_options_iter isa AbstractVector || + block_options_iter isa Dict{Int, Any} + # individual overapproximation options per block + overapproximate_fun = (i, X) -> overapproximate(X, block_options_iter[i]) + else + # uniform overapproximation options for each block + overapproximate_fun = (i, X) -> overapproximate(X, block_options_iter) + end + push!(args, overapproximate_fun) + + # overapproximate function for inputs + lazy_inputs_interval = options[:lazy_inputs_interval] + if lazy_inputs_interval == lazy_inputs_interval_always + overapproximate_inputs_fun = (k, i, x) -> overapproximate_fun(i, x) + else + # first set in a series + function _f(k, i, x::LinearMap{NUM}) where {NUM} + @assert k == 1 "a LinearMap is only expected in the first iteration" + return CacheMinkowskiSum(LazySet{NUM}[x]) + end + # further sets of the series + function _f(k, i, x::MinkowskiSum{NUM, <:CacheMinkowskiSum}) where NUM + if has_constant_directions(block_options_iter, i) + # forget sets if we do not use epsilon-close approximation + forget_sets!(x.X) + end + push!(array(x.X), x.Y) + if lazy_inputs_interval(k) + # overapproximate lazy set + y = overapproximate_fun(i, x.X) + return CacheMinkowskiSum(LazySet{NUM}[y]) + end + return x.X + end + function _f(k, i, x) + # other set types + if lazy_inputs_interval(k) + # overapproximate lazy set + return overapproximate_fun(i, x.X) + end + return x + end + overapproximate_inputs_fun = _f + end + push!(args, overapproximate_inputs_fun) + + # ambient dimension + push!(args, n) + + # number of computed sets + push!(args, N) + + # output function + push!(args, output_function) + + # add mode-specific block(s) argument + push!(args, blocks) + push!(args, partition) + + # time step + push!(args, options[:δ]) + + # termination function + invariant = stateset(problem.s) + vars = [e for block in blocks for e in partition[block]] + if invariant isa Universe + invariant = Universe(length(vars)) + else + invariant = LazySets.Approximations.project(invariant, vars) + end + termination = get_termination_function_out(N, invariant) + push!(args, options[:guards_proj]) + push!(args, options[:block_options_init]) + push!(args, vars) + push!(args, termination) + + info("- Computing successors") + + # progress meter + progress_meter = (N != nothing) ? + Progress(N, 1, "Computing successors ") : + nothing + push!(args, progress_meter) + + # preallocated result + push!(args, res) + + # call the adequate function with the given arguments list + @timing begin + index, skip = reach_blocks!(args...) + end + + # shrink result array + if skip || index < N + if N != nothing + info("termination after only $index of $N steps") + end + deleteat!(res, (skip ? index : index + 1):length(res)) + end + + # return the result + return res +end + +function reach_mixed(problem::Union{IVP{<:CLCS{NUM}, <:LazySet{NUM}}, + IVP{<:CLCCS{NUM}, <:LazySet{NUM}}}, + options::TwoLayerOptions + )::Vector{<:ReachSet} where {NUM <: Real} + info("Time discretization...") + Δ = @timing discretize(problem, options[:δ], + algorithm=options[:discretization], exp_method=options[:exp_method], + sih_method=options[:sih_method]) + return reach_mixed(Δ, options) +end + +function get_termination_function_out(N::Nothing, invariant::Universe) + termination_function = (k, set, t0) -> termination_unconditioned(set) + warn("no termination condition specified; the reachability analysis will " * + "not terminate") + return termination_function +end + +function get_termination_function_out(N::Int, invariant::Universe) + return (k, set, t0) -> termination_N_out(N, k, t0, set) +end + +function get_termination_function_out(N::Nothing, invariant::LazySet) + return (k, set, t0) -> termination_inv_out(invariant, set, t0) +end + +function get_termination_function_out(N::Int, invariant::LazySet) + return (k, set, t0) -> termination_inv_N_out(N, invariant, k, set, t0) +end + +function termination_unconditioned_out(set) + return (false, false, set) +end + +function termination_N_out(N, k, t0, set) + return (k >= N, false, set) +end + +function termination_inv_out(inv, set, t0) + l_inters = Intersection(set, inv) + if isempty(l_inters) + return (true, true, EmptySet()) + else + return (false, false, l_inters) + end +end + +function termination_inv_N_out(N, inv, k, set, t0) + if k >= N + return (true, false, EmptySet()) + end + l_inters = Intersection(set, inv) + if isempty(l_inters) + return (true, true, EmptySet()) + else + return (false, false, l_inters) + end +end diff --git a/src/ReachSets/ContinuousPost/BFFPS19/reach_blocks.jl b/src/ReachSets/ContinuousPost/BFFPS19/reach_blocks.jl new file mode 100644 index 00000000..cacd085c --- /dev/null +++ b/src/ReachSets/ContinuousPost/BFFPS19/reach_blocks.jl @@ -0,0 +1,301 @@ + + +# sparse +function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, + Xhat0::Vector{<:LazySet{NUM}}, + U::Union{ConstantInput, Nothing}, + overapproximate::Function, + overapproximate_inputs::Function, + n::Int, + N::Union{Int, Nothing}, + output_function::Union{Function, Nothing}, + blocks::AbstractVector{Int}, + partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, + δ::NUM, + guards_proj::Vector{<:LazySet{NUM}}, + block_options, + vars::Vector{Int}, + termination::Function, + progress_meter::Union{Progress, Nothing}, + res::Vector{<:ReachSet} + )::Tuple{Int, Bool} where {NUM} + update!(progress_meter, 1) + array = CartesianProductArray(Xhat0[blocks]) + X_store = (output_function == nothing) ? + array : + box_approximation(output_function(array)) + t0 = zero(δ) + t1 = δ + + diff_blocks = setdiff(collect(1:length(partition)),blocks) + + b = length(blocks) + bd = length(diff_blocks) + + array_d = CartesianProductArray(Xhat0[diff_blocks]) + X_store_d = (output_function == nothing) ? + array_d : + box_approximation(output_function(array_d)) + X_store = combine_cpas(X_store, X_store_d, blocks, diff_blocks) + store!(res, 1, X_store, t0, t1, N) + terminate, skip = termination(1, X_store, t0) + if terminate + return 1, skip + end + + b = length(blocks) + Xhatk = Vector{LazySet{NUM}}(undef, b) + Xhatk_d = Vector{LazySet{NUM}}(undef, bd) + ϕpowerk = copy(ϕ) + + if U != nothing + Whatk_blocks = Vector{LazySet{NUM}}(undef, b) + inputs = next_set(U) + @inbounds for i in 1:b + bi = partition[blocks[i]] + Whatk_blocks[i] = overapproximate_inputs(1, blocks[i], proj(bi, n) * inputs) + end + end + + if U != nothing + Whatk_diff_blocks = Vector{LazySet{NUM}}(undef, bd) + inputs = next_set(U) + @inbounds for i in 1:bd + bi = partition[diff_blocks[i]] + Whatk_diff_blocks[i] = overapproximate_inputs(1, diff_blocks[i], proj(bi, n) * inputs) + end + end + + k = 2 + @inbounds while true + update!(progress_meter, k) + for i in 1:b + bi = partition[blocks[i]] + Xhatk_bi = ZeroSet(length(bi)) + for (j, bj) in enumerate(partition) + block = ϕpowerk[bi, bj] + if !iszero(block) + Xhatk_bi = Xhatk_bi + block * Xhat0[j] + end + end + Xhatk_bi_lazy = (U == nothing ? Xhatk_bi : Xhatk_bi + Whatk_blocks[i]) + Xhatk[i] = (output_function == nothing) ? + overapproximate(blocks[i], Xhatk_bi_lazy) : + Xhatk_bi_lazy + end + array = CartesianProductArray(copy(Xhatk)) + X_store = (output_function == nothing) ? + array : + box_approximation(output_function(array)) + + + terminate, skip, rs = termination(k, X_store, t0) + + if !(isdisjoint(X_store, UnionSetArray(guards_proj))) + for i in 1:bd + bi = partition[diff_blocks[i]] + Xhatk_bi = ZeroSet(length(bi)) + for (j, bj) in enumerate(partition) + block = ϕpowerk[bi, bj] + if !iszero(block) + Xhatk_bi = Xhatk_bi + block * Xhat0[j] + end + end + Xhatk_bi_lazy = (U == nothing ? Xhatk_bi : Xhatk_bi + Whatk_diff_blocks[i]) + Xhatk[i] = (output_function == nothing) ? + overapproximate(diff_blocks[i], Xhatk_bi_lazy) : + Xhatk_bi_lazy + end + + array_d = CartesianProductArray(copy(Xhatk_d)) + + X_store_d = (output_function == nothing) ? + array_d : + box_approximation(output_function(array_d)) + if !terminate + X_store = combine_cpas(LazySets.Approximations.overapproximate(rs, CartesianProductArray, block_options), X_store_d, blocks, diff_blocks) + end + end + + t0 = t1 + t1 += δ + store!(res, k, X_store, t0, t1, N) + + terminate, skip = termination(k, X_store, t0) + if terminate + break + end + + if U != nothing + for i in 1:b + bi = partition[blocks[i]] + Whatk_blocks[i] = overapproximate_inputs(k, blocks[i], + Whatk_blocks[i] + row(ϕpowerk, bi) * inputs) + end + end + + ##for diffs + if U != nothing + for i in 1:bd + bi = partition[diff_blocks[i]] + Whatk_diff_blocks[i] = overapproximate_inputs(k, diff_blocks[i], + Whatk_diff_blocks[i] + row(ϕpowerk, bi) * inputs) + end + end + + ϕpowerk = ϕpowerk * ϕ + k += 1 + end + + return k, skip +end + +# dense +function reach_blocks!(ϕ::AbstractMatrix{NUM}, + Xhat0::Vector{<:LazySet{NUM}}, + U::Union{ConstantInput, Nothing}, + overapproximate::Function, + overapproximate_inputs::Function, + n::Int, + N::Union{Int, Nothing}, + output_function::Union{Function, Nothing}, + blocks::AbstractVector{Int}, + partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, + δ::NUM, + guards_proj::Vector{<:LazySet{NUM}}, + block_options, + vars::Vector{Int}, + termination::Function, + progress_meter::Union{Progress, Nothing}, + res::Vector{<:ReachSet} + )::Tuple{Int, Bool} where {NUM} + update!(progress_meter, 1) + array = CartesianProductArray(Xhat0[blocks]) + X_store = (output_function == nothing) ? + array : + box_approximation(output_function(array)) + t0 = zero(δ) + t1 = δ + + diff_blocks = setdiff(collect(1:length(partition)),blocks) + + b = length(blocks) + bd = length(diff_blocks) + + array_d = CartesianProductArray(Xhat0[diff_blocks]) + X_store_d = (output_function == nothing) ? + array_d : + box_approximation(output_function(array_d)) + + terminate, skip = termination(1, X_store, t0) + X_store = combine_cpas(X_store, X_store_d, blocks, diff_blocks) + store!(res, 1, X_store, t0, t1, N) + if terminate + return 1, skip + end + + Xhatk = Vector{LazySet{NUM}}(undef, b) + Xhatk_d = Vector{LazySet{NUM}}(undef, bd) + ϕpowerk = copy(ϕ) + ϕpowerk_cache = similar(ϕ) + + if U != nothing + Whatk_blocks = Vector{LazySet{NUM}}(undef, b) + inputs = next_set(U) + @inbounds for i in 1:b + bi = partition[blocks[i]] + Whatk_blocks[i] = overapproximate_inputs(1, blocks[i], proj(bi, n) * inputs) + end + end + + if U != nothing + Whatk_diff_blocks = Vector{LazySet{NUM}}(undef, bd) + inputs = next_set(U) + @inbounds for i in 1:bd + bi = partition[diff_blocks[i]] + Whatk_diff_blocks[i] = overapproximate_inputs(1, diff_blocks[i], proj(bi, n) * inputs) + end + end + + arr_length = (U == nothing) ? length(partition) : length(partition) + 1 + arr = Vector{LazySet{NUM}}(undef, arr_length) + k = 2 + @inbounds while true + update!(progress_meter, k) + for i in 1:b + bi = partition[blocks[i]] + for (j, bj) in enumerate(partition) + arr[j] = ϕpowerk[bi, bj] * Xhat0[j] + end + if U != nothing + arr[arr_length] = Whatk_blocks[i] + end + Xhatk[i] = (output_function == nothing) ? + overapproximate(blocks[i], MinkowskiSumArray(arr)) : + MinkowskiSumArray(copy(arr)) + end + array = CartesianProductArray(copy(Xhatk)) + + X_store = (output_function == nothing) ? + array : + box_approximation(output_function(array)) + + terminate, skip, rs = termination(k, X_store, t0) + if !(isdisjoint(X_store, UnionSetArray(guards_proj))) + for i in 1:bd + bi = partition[diff_blocks[i]] + for (j, bj) in enumerate(partition) + arr[j] = ϕpowerk[bi, bj] * Xhat0[j] + end + if U != nothing + arr[arr_length] = Whatk_diff_blocks[i] + end + Xhatk_d[i] = (output_function == nothing) ? + overapproximate(diff_blocks[i], MinkowskiSumArray(arr)) : + MinkowskiSumArray(copy(arr)) + end + + array_d = CartesianProductArray(copy(Xhatk_d)) + + X_store_d = (output_function == nothing) ? + array_d : + box_approximation(output_function(array_d)) + if !terminate + rs_oa = LazySets.Approximations.overapproximate(rs, CartesianProductArray, block_options) + X_store = combine_cpas(rs_oa, X_store_d, blocks, diff_blocks) + end + end + + t0 = t1 + t1 += δ + store!(res, k, X_store, t0, t1, N) + + if terminate + break + end + + if U != nothing + for i in 1:b + bi = partition[blocks[i]] + Whatk_blocks[i] = overapproximate_inputs(k, blocks[i], + Whatk_blocks[i] + row(ϕpowerk, bi) * inputs) + end + end + + ##for diffs + if U != nothing + for i in 1:bd + bi = partition[diff_blocks[i]] + Whatk_diff_blocks[i] = overapproximate_inputs(k, diff_blocks[i], + Whatk_diff_blocks[i] + row(ϕpowerk, bi) * inputs) + end + end + + _A_mul_B!(ϕpowerk_cache, ϕpowerk, ϕ) + copyto!(ϕpowerk, ϕpowerk_cache) + + k += 1 + end + + return k, skip +end diff --git a/src/ReachSets/DiscretePost/DecomposedDiscretePost.jl b/src/ReachSets/DiscretePost/DecomposedDiscretePost.jl new file mode 100644 index 00000000..7909f0ad --- /dev/null +++ b/src/ReachSets/DiscretePost/DecomposedDiscretePost.jl @@ -0,0 +1,144 @@ +export DecomposedDiscretePost + +""" + DecomposedDiscretePost <: DiscretePost + +Textbook implementation of a discrete post operator, but with lazy decomposed intersections. + +### Fields + +- `options` -- an `Options` structure that holds the algorithm-specific options + +### Algorithm + +The algorithm is based on [Flowpipe-Guard Intersection for Reachability +Computations with Support Functions](http://spaceex.imag.fr/sites/default/files/frehser_adhs2012.pdf). +""" +struct DecomposedDiscretePost <: DiscretePost + options::Options + + function DecomposedDiscretePost(𝑂::Options) + 𝑂copy = copy(𝑂) + # TODO: Check why it takes always default value for convex_hull + check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:overapproximation], Hyperrectangle) + check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:out_vars], Vector{Int}()) + + return new(𝑂copy) + end +end + +# convenience constructor from pairs of symbols +DecomposedDiscretePost(𝑂::Pair{Symbol,<:Any}...) = DecomposedDiscretePost(Options(Dict{Symbol,Any}(𝑂))) + +# default options for the DecomposedDiscretePost discrete post operator +DecomposedDiscretePost() = DecomposedDiscretePost(Options()) + +init(𝒫::DecomposedDiscretePost, 𝒮::AbstractSystem, 𝑂::Options) = init!(𝒫, 𝒮, copy(𝑂)) + +# TODO: use 𝑂 only? +function init!(𝒫::DecomposedDiscretePost, 𝒮::AbstractSystem, 𝑂::Options) + 𝑂[:n] = statedim(𝒮, 1) + + # solver-specific options (adds default values for unspecified options) + 𝑂out = validate_solver_options_and_add_default_values!(𝑂) + + return 𝑂out +end + +function tube⋂inv!(𝒫::DecomposedDiscretePost, + reach_tube::Vector{<:ReachSet{<:LazySet, N}}, + invariant, + Rsets, + start_interval + ) where {N} + + dirs = 𝒫.options[:overapproximation] + + # counts the number of sets R⋂I added to Rsets + count = 0 + @inbounds for reach_set in reach_tube + push!(Rsets, ReachSet{LazySet{N}, N}(reach_set.X, + reach_set.t_start + start_interval[1], + reach_set.t_end + start_interval[2])) + count = count + 1 + end + + return count +end + +function post(𝒫::DecomposedDiscretePost, + HS::HybridSystem, + waiting_list::Vector{Tuple{Int, ReachSet{LazySet{N}, N}, Int}}, + passed_list, + source_loc_id, + tube⋂inv, + count_Rsets, + jumps, + options + ) where {N} + jumps += 1 + oa = 𝒫.options[:overapproximation] + temp_vars = 𝒫.options[:temp_vars] + source_invariant = HS.modes[source_loc_id].X + inv_isa_Hrep, inv_isa_H_polytope = get_Hrep_info(source_invariant) + + for trans in out_transitions(HS, source_loc_id) + info("Considering transition: $trans") + target_loc_id = target(HS, trans) + target_loc = HS.modes[target(HS, trans)] + target_invariant = target_loc.X + constrained_map = resetmap(HS, trans) + guard = stateset(constrained_map) + # perform jumps + post_jump = Vector{ReachSet{LazySet{N}, N}}() + sizehint!(post_jump, count_Rsets) + for reach_set in tube⋂inv[length(tube⋂inv) - count_Rsets + 1 : end] + if (dim(reach_set.X) == length(temp_vars)) + continue + end + # check intersection with guard + R⋂G = Intersection(reach_set.X, guard) + if isempty(R⋂G) + continue + end + R⋂G = overapproximate(R⋂G, CartesianProductArray, oa) + + # apply assignment + A⌜R⋂G⌟ = apply_assignment(𝒫, constrained_map, R⋂G) + A⌜R⋂G⌟ = overapproximate(A⌜R⋂G⌟, CartesianProductArray, oa) + + # intersect with target invariant + A⌜R⋂G⌟⋂I = Intersection(target_invariant, A⌜R⋂G⌟) + if isempty(A⌜R⋂G⌟⋂I) + continue + end + + A⌜R⋂G⌟⋂I = overapproximate(A⌜R⋂G⌟⋂I, CartesianProductArray, oa) + + + # store result + push!(post_jump, ReachSet{LazySet{N}, N}(A⌜R⋂G⌟⋂I, + reach_set.t_start, + reach_set.t_end)) + end + + postprocess(𝒫, HS, post_jump, options, waiting_list, passed_list, + target_loc_id, jumps) + end +end + +# --- handling assignments --- + +function apply_assignment(𝒫::DecomposedDiscretePost, + constrained_map::Union{IdentityMap, ConstrainedIdentityMap}, + R⋂G::LazySet; + kwargs...) + return R⋂G +end + +function apply_assignment(𝒫::DecomposedDiscretePost, + constrained_map::ConstrainedLinearMap, + R⋂G::LazySet; + kwargs...) + return LinearMap(constrained_map.A, R⋂G) +end diff --git a/src/ReachSets/ReachSets.jl b/src/ReachSets/ReachSets.jl index 9d912224..8a504a27 100644 --- a/src/ReachSets/ReachSets.jl +++ b/src/ReachSets/ReachSets.jl @@ -102,6 +102,8 @@ include("ContinuousPost/GLGM06/GLGM06.jl") include("ContinuousPost/TMJets/TMJets.jl") +include("ContinuousPost/BFFPS19/BFFPS19.jl") + # ======================== # Reachability Algorithms # ======================== @@ -125,6 +127,7 @@ export available_algorithms # ========================== include("DiscretePost/LazyDiscretePost.jl") include("DiscretePost/ConcreteDiscretePost.jl") +include("DiscretePost/DecomposedDiscretePost.jl") # ============================================== # Projection of the reach set in two dimensions diff --git a/src/solve.jl b/src/solve.jl index 124df149..06a9a93e 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -162,6 +162,14 @@ function solve!(system::InitialValueProblem{<:HybridSystem, # TODO temporary hack, to be resolved in #447 options_copy[:mode] = "reach" end + + if opC isa BFFPS19 + opC.options.specified[:HS] = HS + opC.options.specified[:loc_id] = loc_id + opC.options.specified[:opD] = opD + end + + reach_tube = solve!(IVP(loc, X0.X), options_copy, op=opC) if opC isa BFFPSV18 @@ -201,9 +209,9 @@ function solve!(system::InitialValueProblem{<:HybridSystem, if jumps == max_jumps continue end - post(opD, HS, waiting_list, passed_list, loc_id, Rsets, count_Rsets, - jumps, options) + jumps, options) + end if options[:mode] == "check" return CheckSolution(true, -1, options)