Approximations

Module Approximations.jl – polygonal approximation of convex sets through support vectors.

decompose(S::LazySet{N};
          [set_type]::Type{<:Union{HPolygon, Hyperrectangle, Interval}}=Hyperrectangle,
          [ε]::Real=Inf,
          [blocks]::AbstractVector{Int}=default_block_structure(S, set_type),
          [block_types]::Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}(),
          [directions]::Union{Type{<:AbstractDirections}, Nothing}=nothing
         )::CartesianProductArray where {N<:Real}

Decompose a high-dimensional set into a Cartesian product of overapproximations of the projections over the specified subspaces.

Input

• S – set
• set_type – (optional, default: Hyperrectangle) type of set approximation for each subspace
• ε – (optional, default: Inf) error bound for polytopic approximation
• blocks – (optional, default: [2, …, 2] or [1, …, 1] if set_type is an interval) block structure - a vector with the size of each block
• block_types – (optional, default: Interval for 1D and Hyperrectangle for mD blocks) a mapping from set types to blocks
• directions – (optional, default: nothing) template direction type, or nothing

Output

A CartesianProductArray containing the low-dimensional approximated projections.

Algorithm

For each block a specific project method is called, dispatched on the set_type argument.

Notes

If directions is different from nothing, the template directions are used together with blocks. Otherwise, if block_types is given, the options set_type and blocks are ignored.

Examples

The decompose function supports different options, such as: supplying different dimensions for the decomposition, defining the target set of the decomposition, or specifying the degree of accuracy of the target decomposition. You can also choose to make the approximations in low dimensions using template directions. These options are exemplified below.

Different dimensions

By default, decompose returns a Cartesian product of 2D Hyperrectangle sets. For example:

julia> import LazySets.Approximations:decompose

julia> S = Ball2(zeros(4), 1.);

julia> array(decompose(S))
2-element Array{LazySet{Float64},1}:
 Hyperrectangle{Float64}([0.0, 0.0], [1.0, 1.0])
 Hyperrectangle{Float64}([0.0, 0.0], [1.0, 1.0])

Other block sizes can be specified using the blocks option, which refers to each block size of the partition:

julia> array(decompose(S, blocks=[1, 3]))
2-element Array{LazySet{Float64},1}:
 Hyperrectangle{Float64}([0.0], [1.0])
 Hyperrectangle{Float64}([0.0, 0.0, 0.0], [1.0, 1.0, 1.0])

julia> array(decompose(S, blocks=[4]))
1-element Array{LazySet{Float64},1}:
 Hyperrectangle{Float64}([0.0, 0.0, 0.0, 0.0], [1.0, 1.0, 1.0, 1.0])

Different set types

We can also decompose using polygons in constraint representation, through the set_type optional argument:

julia> all([ai isa HPolygon for ai in array(decompose(S, set_type=HPolygon))])
true

For decomposition into 1D subspaces, we can use Interval:

julia> all([ai isa Interval for ai in array(decompose(S, set_type=Interval))])
true

However, if you need to specify different set types for different blocks, the interface presented so far does not apply. In the paragraph Advanced different set types input we explain the input block_types, that can be used precisely for that purpose.

Refining the decomposition I: $ε$-close approximation

The $ε$ option can be used to refine, that is obtain a more accurate decomposition in those blocks where HPolygon types are used, and it relies on the iterative refinement algorithm provided in the Approximations module.

To illustrate this, consider the unit 4D ball in the 2-norm. Using smaller $ε$ implies a better precision, thus more constraints in each 2D decomposition:

julia> S = Ball2(zeros(4), 1.);

julia> d(ε, bi) = array(decompose(S, set_type=HPolygon, ε=ε))[bi]
d (generic function with 1 method)

julia> [length(constraints_list(d(ε, 1))) for ε in [Inf, 0.1, 0.01]]
3-element Array{Int64,1}:
  4
  8
 32

Refining the decomposition II: template polyhedra

Another way to refine the decomposition is using template polyhedra. The idea is to specify a set of template directions, and on each block, compute the polytopic overapproximation obtained by evaluating the support function of the given input set over the template directions.

For example, octagonal 2D approximations of the ball S are obtained with:

julia> B = decompose(S, directions=OctDirections);

julia> length(B.array) == 2 && all(dim(bi) == 2 for bi in B.array)
true

See template_directions.jl for the available template directions. Note that, in contrast to the polygonal $ε$-close approximation, this method can be applied for blocks of any size.

julia> B = decompose(S, directions=OctDirections, blocks=[4]);

julia> length(B.array) == 1 && dim(B.array[1]) == 4
true

Advanced different set types input

We can define different set types for different blocks, using the optional block_types input argument. It is a dictionary where the keys correspond to set types, and the values correspond to the blocks, namely the initial and final block indices should be given.

For example:

julia> S = Ball2(zeros(3), 1.);

julia> array(decompose(S, block_types=Dict(Interval=>[1:1], Hyperrectangle=>[2:3])))
2-element Array{LazySet{Float64},1}:
 Interval{Float64,IntervalArithmetic.Interval{Float64}}([-1, 1])
 Hyperrectangle{Float64}([0.0, 0.0], [1.0, 1.0])

We can additionally pass ε, which is automatically used for each HPolygon type block.

julia> S = Ball2(zeros(8), 1.);

julia> bt = Dict(Interval=>[1:1], Hyperrectangle=>[2:4], HPolygon=>[5:6, 7:8]);

julia> [typeof(ai) for ai in array(decompose(S, block_types=bt, ε=0.01))]
4-element Array{DataType,1}:
 Interval{Float64,IntervalArithmetic.Interval{Float64}}
 Hyperrectangle{Float64}
 HPolygon{Float64}
 HPolygon{Float64}

default_block_structure(S::LazySet, set_type::Type{<:LazySet})::AbstractVector{Int}

Compute the default block structure.

Input

• S – set
• set_type – target set type

Output

A vector representing the block structure, such that:

• If the target set_type is an interval, the default is blocks of size 1.
• Otherwise, the default is blocks of size 2. Depending on the dimension, the last block has size 1 or 2.

project(S::LazySet{N},
        block::AbstractVector{Int},
        set_type::Type{<:LazySet},
        [n]::Int=dim(S),
        [ε]::Real=Inf
       )::LazySet{N} where {N<:Real}

Default implementation for projecting a high-dimensional set to a given set type with possible overapproximation.

Input

• S – set
• block – block structure - a vector with the dimensions of interest
• set_type – target set type
• n – (optional, default: dim(S)) ambient dimension of the set S
• ε – (optional, default: Inf) ignored

Output

A set of type set_type representing an overapproximation of the projection of S.

Algorithm

1. Project the set S with M⋅S, where M is the identity matrix in the block coordinates and zero otherwise.
2. Overapproximate the projected lazy set using overapproximate.

project(S::LazySet{N},
        block::AbstractVector{Int},
        set_type::Type{<:HPolygon},
        [n]::Int=dim(S),
        [ε]::Real=Inf
       )::HPolygon where {N<:Real}

Project a high-dimensional set to a two-dimensional polygon with a certified error bound.

Input

• S – set
• block – block structure - a vector with the dimensions of interest
• set_type – HPolygon - used for dispatch
• n – (optional, default: dim(S)) ambient dimension of the set S
• ε – (optional, default: Inf) error bound for polytopic approximation

Output

A HPolygon representing the epsilon-close approximation of the box approximation of the projection of S.

Notes

block must have length 2.

Algorithm

If ε < Inf, the algorithm proceeds as follows:

1. Project the set S with M⋅S, where M is the identity matrix in the block coordinates and zero otherwise.
2. Overapproximate the set with the given error bound ε.

If ε == Inf, the algorithm uses a box approximation.

project(S::LazySet{N},
        block::AbstractVector{Int},
        set_type::Type{<:Hyperrectangle},
        [n]::Int=dim(S),
        [ε]::Real=Inf
       )::Hyperrectangle where {N<:Real}

Project a high-dimensional set to a low-dimensional hyperrectangle.

Input

• S – set
• block – block structure - a vector with the dimensions of interest
• set_type – Hyperrectangle - used for dispatch
• n – (optional, default: dim(S)) ambient dimension of the set S
• ε – (optional, default: Inf) - used for dispatch, ignored

Output

The box approximation of the projection of S.

project(S::LazySet{N},
        block::AbstractVector{Int},
        directions::AbstractDirections{N},
        n::Int
       )::HPolytope where {N<:Real}

Project a high-dimensional set to a low-dimensional set using template directions.

Input

• S – set
• block – block structure - a vector with the dimensions of interest
• directions – template directions
• n – (optional, default: dim(S)) ambient dimension of the set S

Output

The template direction approximation of the projection of S.

function project(S::LazySet{N},
                block::AbstractVector{Int},
                set_type::Type{<:LinearMap},
                n::Int=dim(S),
                ε::Real=Inf)::LinearMap

Project a high-dimensional set to a low-dimensional set by a lazy linear map.

Input

• S – set
• block – block structure - a vector with the dimensions of interest
• set_type – Hyperrectangle - used for dispatch
• n – (optional, default: dim(S)) ambient dimension of the set S
• ε – (optional, default: Inf) - used for dispatch, ignored

Output

A lazy LinearMap representing a projection of the high-dimensional set to a low-dimensional.

overapproximate(X::S, ::Type{S}) where {S<:LazySet}

Overapproximating a set of type S with type S is a no-op.

Input

• X – set
• Type{S} – set type

Output

The input set.

overapproximate(S::LazySet{N},
                ::Type{<:HPolygon},
                [ε]::Real=Inf)::HPolygon where {N<:Real}

Return an approximation of a given 2D convex set. If no error tolerance is given, or is Inf, the result is a box-shaped polygon. Otherwise the result is an ε-close approximation as a polygon.

Input

• S – convex set, assumed to be two-dimensional
• HPolygon – type for dispatch
• ε – (optional, default: Inf) error bound

Output

A polygon in constraint representation.

overapproximate(S::LazySet, ε::Real)::HPolygon

Alias for overapproximate(S, HPolygon, ε).

overapproximate(S::LazySet,
                Type{<:Hyperrectangle})::Union{Hyperrectangle, EmptySet}

Return an approximation of a given set as a hyperrectangle.

Input

• S – set
• Hyperrectangle – type for dispatch

Output

A hyperrectangle.

overapproximate(S::LazySet)::Union{Hyperrectangle, EmptySet}

Alias for overapproximate(S, Hyperrectangle).

overapproximate(S::ConvexHull{N, Zonotope{N}, Zonotope{N}},
                ::Type{<:Zonotope})::Zonotope where {N<:Real}

Overapproximate the convex hull of two zonotopes.

Input

• S – convex hull of two zonotopes of the same order
• Zonotope – type for dispatch

Algorithm

This function implements the method proposed in Reachability of Uncertain Linear Systems Using Zonotopes, A. Girard, HSCC 2005. The convex hull of two zonotopes of the same order, that we write $Z_j = ⟨c^{(j)}, g^{(j)}_1, …, g^{(j)}_p⟩$ for $j = 1, 2$, can be overapproximated as follows:

It should be noted that the output zonotope is not necessarily the minimal enclosing zonotope, which is in general expensive in high dimensions. This is further investigated in: Zonotopes as bounding volumes, L. J. Guibas et al, Proc. of Symposium on Discrete Algorithms, pp. 803-812.

overapproximate(Z::Zonotope, ::Type{<:Hyperrectangle})::Hyperrectangle

Return a tight overapproximation of a zonotope with an axis-aligned box.

Input

• Z – zonotope
• Hyperrectangle – type for dispatch

Output

A hyperrectangle.

Algorithm

This function implements the method in [Section 5.1.2, 1]. A zonotope $Z = ⟨c, G⟩$ can be overapproximated tightly by an axis-aligned box (i.e. a Hyperrectangle) such that its center is $c$ and the radius along dimension $i$ is the column-sum of the absolute values of the $i$-th row of $G$ for $i = 1,…, p$, where $p$ is the number of generators of $Z$.

[1] Althoff, M., Stursberg, O., & Buss, M. (2010). Computing reachable sets of hybrid systems using a combination of zonotopes and polytopes. Nonlinear analysis: hybrid systems, 4(2), 233-249.

overapproximate(X::LazySet{N}, dir::AbstractDirections{N})::HPolytope{N}
    where {N}

Overapproximating a set with template directions.

Input

• X – set
• dir – direction representation

Output

A HPolytope overapproximating the set X with the directions from dir.

overapproximate(S::LazySet{N}, ::Type{Interval}) where {N<:Real}

Return the overapproximation of a real unidimensional set with an interval.

Input

• S – one-dimensional set
• Interval – type for dispatch

Output

An interval.

overapproximate(cap::Intersection{N,
                                  <:LazySet,
                                  <:Union{AbstractPolytope{N}, HPolyhedron{N}}},
                dir::AbstractDirections{N};
                kwargs...
               ) where {N<:Real}

Return the overapproximation of the intersection between a compact set and a polytope given a set of template directions.

Input

• cap – intersection of a compact set and a polytope
• dir – template directions
• kwargs – additional arguments that are passed to the support function algorithm

Output

A polytope in H-representation such that the normal direction of each half-space is given by an element of dir.

Algorithm

Let di be a direction drawn from the set of template directions dir. Let X be the compact set and let P be the polytope. We overapproximate the set X ∩ H with a polytope in constraint representation using a given set of template directions dir.

The idea is to solve the univariate optimization problem ρ(di, X ∩ Hi) for each half-space in the set P and then take the minimum. This gives an overapproximation of the exact support function.

This algorithm is inspired from G. Frehse, R. Ray. Flowpipe-Guard Intersection for Reachability Computations with Support Functions.

Notes

This method relies on having available the constraints_list of the polytope P.

This method of overapproximations can return a non-empty set even if the original intersection is empty.

overapproximate(cap::Intersection{N, <:HalfSpace{N}, <:AbstractPolytope{N}},
                dir::AbstractDirections{N};
                [kwargs]...
               ) where {N<:Real}

Return the overapproximation of the intersection between a half-space and a polytope given a set of template directions.

Input

• cap – intersection of a half-space and a polytope
• dir – template directions
• kwargs – additional arguments that are passed to the support function algorithm

Output

A polytope in H-representation such that the normal direction of each half-space is given by an element of dir.

ballinf_approximation(S::LazySet{N};
                     )::BallInf{N} where {N<:Real}

Overapproximate a convex set by a tight ball in the infinity norm.

Input

• S – convex set

Output

A tight ball in the infinity norm.

Algorithm

The center and radius of the box are obtained by evaluating the support function of the given convex set along the canonical directions.

box_approximation(S::LazySet{N})::Union{Hyperrectangle{N}, EmptySet{N}}
    where {N<:Real}

Overapproximate a convex set by a tight hyperrectangle.

Input

• S – convex set

Output

A tight hyperrectangle.

Algorithm

The center of the hyperrectangle is obtained by averaging the support function of the given set in the canonical directions, and the lengths of the sides can be recovered from the distance among support functions in the same directions.

interval_hull

Alias for box_approximation.

box_approximation_symmetric(S::LazySet{N}
                           )::Union{Hyperrectangle{N}, EmptySet{N}}
                            where {N<:Real}

Overapproximate a convex set by a tight hyperrectangle centered in the origin.

Input

• S – convex set

Output

A tight hyperrectangle centered in the origin.

Algorithm

The center of the box is the origin, and the radius is obtained by computing the maximum value of the support function evaluated at the canonical directions.

symmetric_interval_hull

Alias for box_approximation_symmetric.

box_approximation_helper(S::LazySet{N};
                        ) where {N<:Real}

Common code of box_approximation and box_approximation_symmetric.

Input

• S – convex set

Output

A tuple containing the data that is needed to construct a tightly overapproximating hyperrectangle.

• c – center
• r – radius

Algorithm

The center of the hyperrectangle is obtained by averaging the support function of the given convex set in the canonical directions. The lengths of the sides can be recovered from the distance among support functions in the same directions.

LocalApproximation{N<:Real}

Type that represents a local approximation in 2D.

Fields

• p1 – first inner point
• d1 – first direction
• p2 – second inner point
• d2 – second direction
• q – intersection of the lines l1 ⟂ d1 at p1 and l2 ⟂ d2 at p2
• refinable – states if this approximation is refinable
• err – error upper bound

Notes

The criteria for being refinable are determined in the method new_approx.

PolygonalOverapproximation{N<:Real}

Type that represents the polygonal approximation of a convex set.

Fields

• S – convex set
• approx_stack – stack of local approximations that still need to be examined
• constraints – vector of linear constraints that are already finalized (i.e., they satisfy the given error bound)

new_approx(S::LazySet, p1::Vector{N}, d1::Vector{N}, p2::Vector{N},
           d2::Vector{N}) where {N<:Real}

Create a LocalApproximation instance for the given excerpt of a polygonal approximation.

Input

• S – convex set
• p1 – first inner point
• d1 – first direction
• p2 – second inner point
• d2 – second direction

Output

A local approximation of S in the given directions.

addapproximation!(Ω::PolygonalOverapproximation, p1::Vector{N},
    d1::Vector{N}, p2::Vector{N}, d2::Vector{N}) where {N<:Real}

Input

• Ω – polygonal overapproximation of a convex set
• p1 – first inner point
• d1 – first direction
• p2 – second inner point
• d2 – second direction

Output

The list of local approximations in Ω of the set Ω.S is updated in-place and the new approximation is returned by this function.

refine(approx::LocalApproximation, S::LazySet
      )::Tuple{LocalApproximation, LocalApproximation}

Refine a given local approximation of the polygonal approximation of a convex set by splitting along the normal direction of the approximation.

Input

• approx – local approximation to be refined
• S – 2D convex set

Output

The tuple consisting of the refined right and left local approximations.

tohrep(Ω::PolygonalOverapproximation{N})::AbstractHPolygon{N}
    where {N<:Real}

Convert a polygonal overapproximation into a concrete polygon.

Input

• Ω – polygonal overapproximation of a convex set

Output

A polygon in constraint representation.

Algorithm

Internally we keep the constraints sorted. Hence we do not need to use addconstraint! when creating the HPolygon.

approximate(S::LazySet{N},
            ε::N)::PolygonalOverapproximation{N} where {N<:Real}

Return an ε-close approximation of the given 2D convex set (in terms of Hausdorff distance) as an inner and an outer approximation composed by sorted local Approximation2D.

Input

• S – 2D convex set
• ε – error bound

Output

An ε-close approximation of the given 2D convex set.

constraint(approx::LocalApproximation)

Convert a local approximation to a linear constraint.

Input

• approx – local approximation

Output

A linear constraint.

AbstractDirections{N}

Abstract type for template direction representations.

Notes

All subtypes should implement the standard iterator methods from Base and the function dim(d<:AbstractDirections)::Int.

BoxDirections{N} <: AbstractDirections{N}

Box direction representation.

Fields

• n – dimension

OctDirections{N} <: AbstractDirections{N}

Octagon direction representation.

Fields

• n – dimension

Notes

Octagon directions consist of all vectors that are zero almost everywhere except in two dimensions $i$, $j$ (possibly $i = j$) where it is $±1$.

BoxDiagDirections{N} <: AbstractDirections{N}

Box-diagonal direction representation.

Fields

• n – dimension

Notes

Box-diagonal directions can be seen as the union of diagonal directions (all entries are ±1) and box directions (one entry is ±1, all other entries are 0). The iterator first enumerates all diagonal directions, and then all box directions.