# Approximations

This section of the manual describes the Cartesian decomposition algorithms and the approximation of high-dimensional convex sets using projections.

`LazySets.Approximations`

— ModuleModule `Approximations.jl`

– polygonal approximation of sets.

## Cartesian Decomposition

`LazySets.Approximations.decompose`

— Function```
decompose(S::LazySet{N},
partition::AbstractVector{<:AbstractVector{Int}},
block_options) where {N}
```

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

**Input**

`S`

– set`partition`

– vector of blocks (i.e., of vectors of integers) (see the Notes below)`block_options`

– mapping from block indices in`partition`

to a corresponding overapproximation option; we only require access via`[⋅]`

(but see also the Notes below)

**Output**

A `CartesianProductArray`

containing the low-dimensional approximated projections.

**Algorithm**

For each block a specific `project`

method is called, dispatching on the corresponding overapproximation option.

**Notes**

The argument `partition`

deserves some explanation. Typically, the list of blocks should form a partition of the set $\{1, \dots, n\}$ represented as a list of consecutive blocks, where $n$ is the ambient dimension of set `S`

.

However, technically there is no problem if the blocks are not consecutive, blocks are missing, blocks occur more than once, or blocks are overlapping. The resulting set must be interpreted with care in such cases (e.g., it will not necessarily be an overapproximation of `S`

).

For convenience, the argument `block_options`

can also be given as a single option instead of a mapping, which is then interpreted as the option for all blocks.

**Examples**

The argument `block_options`

supports different options: one can specify the target set, the degree of accuracy, and template directions. These options are exemplified below, where we use the following example.

```
julia> S = Ball2(zeros(4), 1.0); # set to be decomposed (4D 2-norm unit ball)
julia> P2d = [1:2, 3:4]; # a partition with two blocks, each of size two
julia> P1d = [[1], [2], [3], [4]]; # a partition with four blocks, each of size one
```

**Different set types**

We can decompose using polygons in constraint representation:

```
julia> Y = decompose(S, P2d, HPolygon);
julia> all(ai isa HPolygon for ai in array(Y))
true
```

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

:

```
julia> Y = decompose(S, P1d, Interval);
julia> all(ai isa Interval for ai in array(Y))
true
```

However, if you need to specify different set types for different blocks, the interface presented so far does not apply. See the paragraph *Advanced input for different block approximations* below for how to do that.

**Refining the decomposition: $ε$-close approximation**

The $ε$ option can be used to refine a decomposition, i.e., obtain a more accurate result. We use the Iterative refinement algorithm from the `Approximations`

module.

To illustrate this, consider again the set `S`

from above. We decompose into two 2D polygons. Using smaller $ε$ implies a better precision, thus more constraints in each 2D decomposition. In the following example, we look at the number of constraints in the first block.

```
julia> d(ε, bi) = array(decompose(S, P2d, (HPolygon => ε)))[bi]
d (generic function with 1 method)
julia> [length(constraints_list(d(ε, 1))) for ε in [Inf, 0.1, 0.01]]
3-element Vector{Int64}:
4
8
32
```

**Refining the decomposition: template polyhedra**

Another way to refine a decomposition is by using template polyhedra. The idea is to specify a set of template directions and then compute on each block 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 set `S`

are obtained with:

```
julia> B = decompose(S, P2d, OctDirections);
julia> length(B.array) == 2 && all(dim(bi) == 2 for bi in B.array)
true
```

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

```
julia> B = decompose(S, [1:4], OctDirections);
julia> length(B.array) == 1 && dim(B.array[1]) == 4
true
```

**Advanced input for different block approximations**

Instead of defining the approximation option uniformly for each block, we can define different approximations for different blocks. For this purpose, the argument `block_options`

can also be a mapping from block index (in the partition) to the corresponding approximation option.

For example, we can approximate the first block with a `Hyperrectangle`

and the second block with $ε$-close approximation for $ε = 0.1$:

```
julia> res = array(decompose(S, P2d, Dict(1 => Hyperrectangle, 2 => 0.1)));
julia> typeof(res[1]), typeof(res[2])
(Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}, HPolygon{Float64, Vector{Float64}})
```

`decompose(S::LazySet, block_options; [block_size]::Int=1)`

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

**Input**

`S`

– set`block_options`

– overapproximation option or mapping from block indices to a corresponding overapproximation option`block_size`

– (optional; default:`1`

) size of the blocks

**Output**

A `CartesianProductArray`

containing the low-dimensional approximated projections.

## Overapproximations

`LazySets.Approximations.overapproximate`

— Function`overapproximate(X::S, ::Type{S}, args...) where {S<:LazySet}`

Overapproximating a set of type `S`

with type `S`

is a no-op.

**Input**

`X`

– set`Type{S}`

– target set type`args`

– further arguments (ignored)`kwargs`

– further keyword arguments (ignored)

**Output**

The input set.

`overapproximate(S::LazySet)`

Alias for `overapproximate(S, Hyperrectangle)`

resp. `box_approximation(S)`

.

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

Alias for `box_approximation(S)`

.

`overapproximate(S::LazySet, ::Type{<:BallInf})`

Alias for `ballinf_approximation(S)`

.

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

Overapproximate a given 2D set using iterative refinement.

**Input**

`S`

– two-dimensional bounded set`HPolygon`

– target set type`ε`

– (optional, default:`Inf`

) error tolerance`prune`

– (optional, default:`true`

) flag for removing redundant constraints in the end

**Output**

A polygon in constraint representation.

**Notes**

The result is always a convex overapproximation of the input set.

If no error tolerance ε is given, or is `Inf`

, the result is a box-shaped polygon. For convex input sets, the result is an ε-close polygonal overapproximation with respect to the Hausdorff distance.

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

Alias for `overapproximate(S, HPolygon, ε)`

.

```
overapproximate(X::LazySet{N}, dirs::AbstractDirections;
[prune]::Bool=true) where {N}
```

Overapproximate a (possibly unbounded) set with template directions.

**Input**

`X`

– set`dirs`

– directions`prune`

– (optional, default:`true`

) flag for removing redundant constraints

**Output**

A polyhedron overapproximating the set `X`

with the directions from `dirs`

. The overapproximation is computed using the support function. The result is an `HPolytope`

if it is bounded and otherwise an `HPolyhedron`

.

`overapproximate(X::LazySet{N}, dirs::Type{<:AbstractDirections}) where {N}`

Overapproximate a set with template directions.

**Input**

`X`

– set`dirs`

– type of direction representation

**Output**

A polyhedron overapproximating the set `X`

with the directions from `dirs`

. The result is an `HPolytope`

if it is bounded and otherwise an `HPolyhedron`

.

```
overapproximate(cap::Intersection{N, <:LazySet, <:AbstractPolyhedron},
dirs::AbstractDirections;
kwargs...
) where {N}
```

Overapproximate the intersection between a set and a polyhedron given a set of template directions.

**Input**

`cap`

– intersection of a set and a polyhedron`dirs`

– template directions`kwargs`

– additional arguments that are passed to the support function algorithm

**Output**

A polytope or polyhedron in H-representation such that the normal direction of each half-space is given by an element of `dirs`

.

**Algorithm**

Let `di`

be a direction drawn from the set of template directions `dirs`

. Let `X`

be the set and let `P`

be the polyhedron. We overapproximate the set `X ∩ H`

with a polytope or polyhedron in constraint representation using a given set of template directions `dirs`

.

The idea is to solve the univariate optimization problem `ρ(di, X ∩ Hi)`

for each half-space of the set `P`

and then take the minimum. This gives an overapproximation of the exact support function.

This algorithm is inspired from [1].

**Notes**

This method relies on having available the `constraints_list`

of the polyhedron `P`

.

This method may return a non-empty set even if the original set is empty.

[1] G. Frehse, R. Ray. *Flowpipe-Guard Intersection for Reachability Computations with Support Functions*. ADHS 2012.

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

Overapproximate 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`dirs`

– 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 `dirs`

.

```
overapproximate(Z::AbstractZonotope, ::Type{<:HParallelotope},
[indices]=1:dim(Z))
```

Overapproximate a zonotopic set with a parallelotope in constraint representation.

**Input**

`Z`

– zonotopic set`HParallelotope`

– target set type`indices`

– (optional; default:`1:dim(Z)`

) generator indices selected when constructing the parallelotope

**Output**

An overapproximation of the given zonotopic set using a parallelotope.

**Algorithm**

The algorithm is based on Proposition 8 discussed in Section 5 of [1].

[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::Intersection{N, <:AbstractZonotope, <:Hyperplane},
dirs::AbstractDirections) where {N}
```

Overapproximate the intersection between a zonotopic set and a hyperplane with a polyhedron or polytope using the given directions.

**Input**

`X`

– intersection between a zonotopic set and a hyperplane`dirs`

– type of direction representation

**Output**

An overapproximation of the intersection between a zonotopic set and a hyperplane. If the directions are bounding, the result is an `HPolytope`

, otherwise the result is an `HPolyhedron`

.

**Algorithm**

This function implements [Algorithm 8.1, 1].

[1] Colas Le Guernic. *Reachability Analysis of Hybrid Systems with Linear continuous dynamics* (Doctoral dissertation). 2009.

```
overapproximate(QM::QuadraticMap{N, <:SparsePolynomialZonotope},
::Type{<:SparsePolynomialZonotope}) where {N}
```

Overapproximate a quadratic map of a sparse polynomial zonotope with a sparse polynomial zonotope.

**Input**

`QM`

– quadratic map of a sparse polynomial zonotope`SparsePolynomialZonotope`

– target type

**Output**

A sparse polynomial zonotope overapproximating the quadratic map of a sparse polynomial zonotope.

**Algorithm**

This method implements Proposition 13 of [1].

[1] N. Kochdumper and M. Althoff. *Sparse Polynomial Zonotopes: A Novel Set Representation for Reachability Analysis*. Transactions on Automatic Control

`overapproximate(S::LazySet, ::Type{<:Interval})`

Return the overapproximation of a set with an interval.

**Input**

`S`

– one-dimensional set`Interval`

– target type

**Output**

An interval.

**Algorithm**

We use the `extrema`

function.

`overapproximate(cap::Intersection, ::Type{<:Interval})`

Return the overapproximation of a lazy intersection with an interval.

**Input**

`cap`

– one-dimensional intersection`Interval`

– target type

**Output**

An interval.

**Algorithm**

The algorithm recursively overapproximates the two intersected sets with intervals and then intersects these. (Note that this fails if the sets are unbounded.)

For convex sets this algorithm is precise.

`overapproximate(cap::IntersectionArray, ::Type{<:Interval})`

Return the overapproximation of an intersection array with an interval.

**Input**

`cap`

– one-dimensional intersection array`Interval`

– target type

**Output**

An interval.

**Algorithm**

The algorithm recursively overapproximates the intersected sets with intervals and then intersects these. (Note that this fails if the sets are unbounded.)

For convex sets this algorithm is precise.

`overapproximate(Z::Zonotope, ::Type{<:Zonotope}, r::Union{Integer, Rational})`

Reduce the order of a zonotope.

**Input**

`Z`

– zonotope`Zonotope`

– target set type`r`

– desired order

**Output**

A new zonotope with `r`

generators, if possible.

**Algorithm**

This method falls back to `reduce_order`

with the default algorithm.

```
overapproximate(X::ConvexHull{N, <:AbstractZonotope, <:AbstractZonotope},
::Type{<:Zonotope}) where {N}
```

Overapproximate the convex hull of two zonotopic sets.

**Input**

`X`

– convex hull of two zonotopic sets`Zonotope`

– target set type`algorithm`

– (optional; default:`"mean"`

) choice of algorithm; possible values are`"mean"`

and`"join"`

**Output**

A zonotope $Z$ such that $X ⊆ Z$.

**Algorithm**

The algorithm can be controlled by the parameter `algorithm`

. Note that the results of the two implemented algorithms are generally incomparable.

**'mean' method**

If `algorithm == "mean"`

, we choose the method proposed in [1]. The convex hull of two zonotopic sets $Z₁$ and $Z₂$ of the same order, which we write

\[Z_j = ⟨c^{(j)}, g^{(j)}_1, …, g^{(j)}_p⟩\]

for $j = 1, 2$, can be overapproximated as follows:

\[CH(Z_1, Z_2) ⊆ \frac{1}{2}⟨c^{(1)}+c^{(2)}, g^{(1)}_1+g^{(2)}_1, …, g^{(1)}_p+g^{(2)}_p, c^{(1)}-c^{(2)}, g^{(1)}_1-g^{(2)}_1, …, g^{(1)}_p-g^{(2)}_p⟩.\]

If the zonotope order is not the same, this algorithm calls `reduce_order`

to reduce the order to the minimum of the arguments.

It should be noted that the output zonotope is not necessarily the minimal enclosing zonotope, which is in general expensive to obtain in high dimensions. This is further investigated in [2].

**'join' method**

If `algorithm == "join"`

, we choose the method proposed in [3, Definition 1]. The convex hull $X$ of two zonotopic sets $Z₁$ and $Z₂$ is overapproximated by a zonotope $Z₃$ such that the box approximation of $X$ is identical with the box approximation of $Z₃$. Let $□(X)$ denote the box approximation of $X$. The center of $Z₃$ is the center of $□(X)$.

The generator construction consists of two phases. In the first phase, we construct generators $g$ as a combination of one generator from $Z₁$, say, $g₁$, with another generator from $Z₂$, say, $g₂$. The entry of $g$ in the $i$-th dimension is given as

\[ g[i] = \arg\min_{\min(g₁[i], g₂[i]) ≤ x ≤ \max(g₁[i], g₂[i])} |x|.\]

If $g$ is the zero vector, it can be omitted.

In the second phase, we construct another generator for each dimension. These generators are scaled unit vectors. The following formula defines the sum of all those generators.

\[ \sup(□(X)) - c - ∑_g |g|\]

where $c$ is the center of the new zonotope and the $g$s are the generators constructed in the first phase.

**References**

[1] *Reachability of Uncertain Linear Systems Using Zonotopes*. A. Girard. HSCC 2005.

[2] *Zonotopes as bounding volumes*. L. J. Guibas, A. T. Nguyen, L. Zhang. SODA 2003.

[3] *The zonotope abstract domain Taylor1+*. K. Ghorbal, E. Goubault, S. Putot. CAV 2009.

```
overapproximate(lm::LinearMap{N, <:AbstractZonotope},
::Type{<:Zonotope}) where {N}
```

Overapproximate a lazy linear map of a zonotopic set with a zonotope.

**Input**

`lm`

– lazy linear map of a zonotopic set`Zonotope`

– target set type

**Output**

The tight zonotope corresponding to `lm`

.

```
overapproximate(P::SimpleSparsePolynomialZonotope, ::Type{<:Zonotope};
[nsdiv]=1, [partition]=nothing)
```

Overapproximate a simple sparse polynomial zonotope with a zonotope.

**Input**

`P`

– simple sparse polynomial zonotope`Zonotope`

– target set type`nsdiv`

– (optional, default:`1`

) size of uniform partitioning grid`partition`

– (optional, default:`nothing`

) tuple of integers indicating the number of blocks in each dimensino; the length should match`nparams(P)`

**Output**

A zonotope.

```
overapproximate(P::SimpleSparsePolynomialZonotope, ::Type{<:Zonotope},
dom::IntervalBox)
```

Overapproximate a simple sparse polynomial zonotope over the parameter domain `dom`

with a zonotope.

**Input**

`P`

– simple sparse polynomial zonotope`Zonotope`

– target set type`dom`

– parameter domain, which should be a subset of`[-1, 1]^q`

, where`q = nparams(P)`

**Output**

A zonotope.

`overapproximate(P::SparsePolynomialZonotope, ::Type{<:Zonotope})`

Overapproximate a sparse polynomial zonotope with a zonotope.

**Input**

`P`

– sparse polynomial zonotope`Zonotope`

– target set type

**Output**

A zonotope.

```
overapproximate(X::LazySet, ZT::Type{<:Zonotope},
dir::Union{AbstractDirections, Type{<:AbstractDirections}};
[algorithm]="vrep", kwargs...)
```

Overapproximate a set with a zonotope.

**Input**

`X`

– set`Zonotope`

– target set type`dir`

– directions used for the generators`algorithm`

– (optional, default:`"vrep"`

) algorithm used to compute the overapproximation`kwargs`

– further algorithm-specific options

**Output**

A zonotope that overapproximates `X`

and uses at most the generator directions provided in `dir`

(redundant directions will be ignored).

**Notes**

Two algorithms are available:

`"vrep"`

– Overapproximate a polytopic set with a zonotope of minimal total generator sum using only generators in the given directions. Under this constraint, the zonotope has the minimal sum of generator vectors. See the docstring of`_overapproximate_zonotope_vrep`

for further details.`"cpa"`

– Overapproximate a polytopic set with a zonotope using a Cartesian decomposition into two-dimensional blocks. See the docstring of`_overapproximate_zonotope_cpa`

for further details.

```
overapproximate(r::Rectification{N, <:AbstractZonotope},
::Type{<:Zonotope}) where {N}
```

Overapproximate the rectification of a zonotopic set with a zonotope.

**Input**

`r`

– lazy rectification of a zonotopic set`Zonotope`

– target set type

**Output**

A zonotope overapproximation of the set obtained by rectifying `Z`

.

**Algorithm**

This method implements [Theorem 3.1, 1].

[1] Singh, G., Gehr, T., Mirman, M., Püschel, M., & Vechev, M. *Fast and effective robustness certification*. NeurIPS 2018.

```
overapproximate(CHA::ConvexHullArray{N, <:AbstractZonotope},
::Type{<:Zonotope}) where {N}
```

Overapproximate the convex hull of a list of zonotopic sets with a zonotope.

**Input**

`CHA`

– convex hull array of zonotopic sets`Zonotope`

– target set type

**Output**

A zonotope overapproximation of the convex hull array of zonotopic sets.

**Algorithm**

This method iteratively applies the overapproximation algorithm to the convex hull of two zonotopic sets from the given array of zonotopic sets.

```
overapproximate(QM::QuadraticMap{N, <:AbstractZonotope},
::Type{<:Zonotope}) where {N}
```

Overapproximate a quadratic map of a zonotopic set with a zonotope.

**Input**

`QM`

– quadratic map of a zonotopic set`Zonotope`

– target set type

**Output**

A zonotope overapproximating the quadratic map of a zonotopic set.

**Notes**

Mathematically, a quadratic map of a zonotope with matrices $Q$ is defined as:

\[Z_Q = \right\{ \lambda | \lambda_i = x^T Q\^{(i)} x,~i = 1, \ldots, n,~x \in Z \left\}\]

**Algorithm**

This method implements [Lemma 1, 1].

[1] Matthias Althoff and Bruce H. Krogh. *Avoiding geometric intersection operations in reachability analysis of hybrid systems*. HSCC 2012.

```
overapproximate(X::Intersection{N, <:AbstractZonotope, <:Hyperplane},
::Type{<:Zonotope})
```

Overapproximate the intersection of a zonotopic set and a hyperplane with a zonotope.

**Input**

`X`

– intersection of a zonotopic set and a hyperplane`Zonotope`

– target set type

**Output**

A zonotope overapproximating the intersection.

**Algorithm**

This method implements Algorithm 3 in [1].

[1] Moussa Maïga, Nacim Ramdani, Louise Travé-Massuyès, Christophe Combastel: *A CSP versus a zonotope-based method for solving guard set intersection in nonlinear hybrid reachability*. Mathematics in Computer Science (8) 2014.

```
overapproximate(lm::LinearMap{N, <:CartesianProductArray},
::Type{<:CartesianProductArray{N, S}}
) where {N, S<:LazySet}
```

Decompose a lazy linear map of a Cartesian product array while keeping the original block structure.

**Input**

`lm`

– lazy linear map of Cartesian product array`CartesianProductArray`

– target set type

**Output**

A `CartesianProductArray`

representing the decomposed linear map.

```
overapproximate(lm::LinearMap{N, <:CartesianProductArray},
::Type{<:CartesianProductArray},
dir::Type{<:AbstractDirections}) where {N}
```

Decompose a lazy linear map of a Cartesian product array with template directions while keeping the original block structure.

**Input**

`lm`

– lazy linear map of a Cartesian product array`CartesianProductArray`

– target set type`dir`

– template directions for overapproximation

**Output**

A `CartesianProductArray`

representing the decomposed linear map.

```
overapproximate(lm::LinearMap{N, <:CartesianProductArray},
::Type{<:CartesianProductArray},
set_type::Type{<:LazySet}) where {N}
```

Decompose a lazy linear map of a Cartesian product array with a given set type while keeping the original block structure.

**Input**

`lm`

– lazy linear map of a Cartesian product array`CartesianProductArray`

– target set type`set_type`

– set type for overapproximation

**Output**

A `CartesianProductArray`

representing the decomposed linear map.

```
overapproximate(rm::ResetMap{N, <:CartesianProductArray},
::Type{<:CartesianProductArray}, oa) where {N}
```

Overapproximate a reset map (that only resets to zero) of a Cartesian product with a new Cartesian product.

**Input**

`rm`

– reset map`CartesianProductArray`

– target set type`oa`

– overapproximation option

**Output**

A Cartesian product with the same block structure.

**Notes**

This implementation currently only supports resets to zero.

**Algorithm**

We convert the `ResetMap`

into a `LinearMap`

and then call the corresponding `overapproximate`

method.

```
overapproximate(cap::Intersection{N,
<:CartesianProductArray,
<:AbstractPolyhedron},
::Type{<:CartesianProductArray}, oa) where {N}
```

Overapproximate the intersection of a Cartesian product of a finite number of sets and a polyhedron with a new Cartesian product.

**Input**

`cap`

– lazy intersection of a Cartesian product array and a polyhedron`CartesianProductArray`

– target set type`oa`

– overapproximation option

**Output**

A `CartesianProductArray`

that overapproximates the intersection of `cpa`

and `P`

.

**Algorithm**

The intersection only needs to be computed in the blocks of `cpa`

that are constrained in `P`

. Hence we first collect those constrained blocks in a lower-dimensional Cartesian product array and then convert to an `HPolytope`

`X`

. Then we take the intersection of `X`

and the projection of `Y`

onto the corresponding dimensions. (This projection is purely syntactic and exact.) Finally we decompose the result again and plug together the unaffected old blocks and the newly computed blocks. The result is a `CartesianProductArray`

with the same block structure as in `X`

.

```
overapproximate(vTM::Vector{TaylorModel1{T, S}}, ::Type{<:Zonotope};
[remove_zero_generators]::Bool=true
[normalize]::Bool=true) where {T, S}
```

Overapproximate a Taylor model in one variable with a zonotope.

**Input**

`vTM`

– vector of`TaylorModel1`

`Zonotope`

– target set type`remove_zero_generators`

– (optional; default:`true`

) flag to remove zero generators of the resulting zonotope`normalize`

– (optional; default:`true`

) flag to skip the normalization of the Taylor models

**Output**

A zonotope that overapproximates the range of the given Taylor model.

**Examples**

If the polynomials are linear, this method exactly transforms to a zonotope. The nonlinear case necessarily introduces overapproximation error. Consider the linear case first:

```
julia> using LazySets, TaylorModels
julia> const IA = IntervalArithmetic;
julia> I = IA.Interval(-0.5, 0.5) # interval remainder
[-0.5, 0.5]
julia> x₀ = IA.Interval(0.0) # expansion point
[0, 0]
julia> D = IA.Interval(-3.0, 1.0)
[-3, 1]
julia> p1 = Taylor1([2.0, 1.0], 2) # define a linear polynomial
2.0 + 1.0 t + 𝒪(t³)
julia> p2 = Taylor1([0.9, 3.0], 2) # define another linear polynomial
0.9 + 3.0 t + 𝒪(t³)
julia> vTM = [TaylorModel1(pi, I, x₀, D) for pi in [p1, p2]]
2-element Vector{TaylorModel1{Float64, Float64}}:
2.0 + 1.0 t + [-0.5, 0.5]
0.9 + 3.0 t + [-0.5, 0.5]
```

Here, `vTM`

is a Taylor model vector, since each component is a Taylor model in one variable (`TaylorModel1`

). Using `overapproximate(vTM, Zonotope)`

we can compute its associated zonotope in generator representation:

```
julia> Z = overapproximate(vTM, Zonotope);
julia> center(Z)
2-element Vector{Float64}:
1.0
-2.1
julia> Matrix(genmat(Z))
2×3 Matrix{Float64}:
2.0 0.5 0.0
6.0 0.0 0.5
```

Note how the generators of this zonotope mainly consist of two pieces: one comes from the linear part of the polynomials, and another one corresponds to the interval remainder. This conversion gives the same upper and lower bounds as the range evaluation using interval arithmetic:

```
julia> X = box_approximation(Z)
Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}([1.0, -2.1], [2.5, 6.5])
julia> Y = evaluate(vTM[1], vTM[1].dom) × evaluate(vTM[2], vTM[2].dom)
[-1.5, 3.5] × [-8.60001, 4.40001]
julia> H = convert(Hyperrectangle, Y) # this IntervalBox is the same as X
Hyperrectangle{Float64, StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SVector{2, Float64}}([1.0, -2.1000000000000005], [2.5, 6.500000000000001])
```

However, the zonotope returns better results if we want to approximate the Taylor model because it is not axis-aligned:

```
julia> d = [-0.35, 0.93];
julia> ρ(d, Z) < ρ(d, X)
true
```

This method also works if the polynomials are non-linear; for example suppose that we add a third polynomial with a quadratic term:

```
julia> p3 = Taylor1([0.9, 3.0, 1.0], 3)
0.9 + 3.0 t + 1.0 t² + 𝒪(t⁴)
julia> vTM = [TaylorModel1(pi, I, x₀, D) for pi in [p1, p2, p3]]
3-element Vector{TaylorModel1{Float64, Float64}}:
2.0 + 1.0 t + [-0.5, 0.5]
0.9 + 3.0 t + [-0.5, 0.5]
0.9 + 3.0 t + 1.0 t² + [-0.5, 0.5]
julia> Z = overapproximate(vTM, Zonotope);
julia> center(Z)
3-element Vector{Float64}:
1.0
-2.1
2.4
julia> Matrix(genmat(Z))
3×4 Matrix{Float64}:
2.0 0.5 0.0 0.0
6.0 0.0 0.5 0.0
6.0 0.0 0.0 5.0
```

The last generator corresponds to the addition of the interval remainder and the box overapproximation of the nonlinear part of `p3`

over the domain.

**Algorithm**

Let $\text{vTM} = (p, I)$ be a vector of $m$ Taylor models, where $I$ is the interval remainder in $\mathbb{R}^m$. Let $p_{lin}$ (resp. $p_{nonlin}$) correspond to the linear (resp. nonlinear) part of each scalar polynomial.

The range of $\text{vTM}$ can be enclosed by a zonotope with center $c$ and matrix of generators $G$, $Z = ⟨c, G⟩$, by performing a conservative linearization of $\text{vTM}$:

\[ vTM' = (p', I') := (p_{lin} − p_{nonlin} , I + \text{Int}(p_{nonlin})).\]

This algorithm proceeds in two steps:

1- Conservatively linearize $\text{vTM}$ as above and compute a box overapproximation of the nonlinear part. 2- Transform the linear Taylor model to a zonotope exactly through variable normalization onto the symmetric intervals $[-1, 1]$.

```
overapproximate(vTM::Vector{TaylorModelN{N, T, S}}, ::Type{<:Zonotope};
[remove_zero_generators]::Bool=true
[normalize]::Bool=true) where {N, T, S}
```

Overapproximate a multivariate Taylor model with a zonotope.

**Input**

`vTM`

– vector of`TaylorModelN`

`Zonotope`

– target set type`remove_zero_generators`

– (optional; default:`true`

) flag to remove zero generators of the resulting zonotope`normalize`

– (optional; default:`true`

) flag to skip the normalization of the Taylor models

**Output**

A zonotope that overapproximates the range of the given Taylor model.

**Examples**

Consider a vector of two 2-dimensional Taylor models of order 2 and 4, respectively.

```
julia> using LazySets, TaylorModels
julia> const IA = IntervalArithmetic;
julia> x₁, x₂ = set_variables(Float64, ["x₁", "x₂"], order=8)
2-element Vector{TaylorN{Float64}}:
1.0 x₁ + 𝒪(‖x‖⁹)
1.0 x₂ + 𝒪(‖x‖⁹)
julia> x₀ = IntervalBox(0..0, 2) # expansion point
[0, 0]²
julia> Dx₁ = IA.Interval(0.0, 3.0) # domain for x₁
[0, 3]
julia> Dx₂ = IA.Interval(-1.0, 1.0) # domain for x₂
[-1, 1]
julia> D = Dx₁ × Dx₂ # take the Cartesian product of the domain on each variable
[0, 3] × [-1, 1]
julia> r = IA.Interval(-0.5, 0.5) # interval remainder
[-0.5, 0.5]
julia> p1 = 1 + x₁^2 - x₂
1.0 - 1.0 x₂ + 1.0 x₁² + 𝒪(‖x‖⁹)
julia> p2 = x₂^3 + 3x₁^4 + x₁ + 1
1.0 + 1.0 x₁ + 1.0 x₂³ + 3.0 x₁⁴ + 𝒪(‖x‖⁹)
julia> vTM = [TaylorModelN(pi, r, x₀, D) for pi in [p1, p2]]
2-element Vector{TaylorModelN{2, Float64, Float64}}:
1.0 - 1.0 x₂ + 1.0 x₁² + [-0.5, 0.5]
1.0 + 1.0 x₁ + 1.0 x₂³ + 3.0 x₁⁴ + [-0.5, 0.5]
julia> Z = overapproximate(vTM, Zonotope);
julia> center(Z)
2-element Vector{Float64}:
5.5
124.0
julia> Matrix(genmat(Z))
2×4 Matrix{Float64}:
0.0 -1.0 5.0 0.0
1.5 0.0 0.0 123.0
```

**Algorithm**

We refer to the algorithm description for the univariate case.

`LazySets.Approximations._overapproximate_zonotope_vrep`

— Function```
_overapproximate_zonotope_vrep(X::LazySet{N},
dir::AbstractDirections;
solver=default_lp_solver(N)) where {N}
```

Overapproximate a polytopic set with a zonotope of minimal total generator sum using only generators in the given directions.

**Input**

`X`

– polytopic set`dir`

– directions used for the generators`solver`

– (optional, default:`default_lp_solver(N)`

) the backend used to solve the linear program

**Output**

A zonotope that overapproximates `X`

and uses at most the directions provided in `dir`

(redundant directions will be ignored). Under this constraint, the zonotope has the minimal sum of generator vectors.

**Notes**

The algorithm only requires one representative of each generator direction and their additive inverse (e.g., only one of `[1, 0]`

and `[-1, 0]`

) and assumes that the directions are normalized. We preprocess the directions in that respect.

**Algorithm**

We solve a linear program parametric in the vertices $v_j$ of `X`

and the directions $d_k$ in `dir`

presented in Section 4.2 in [1], adapting the notation to the one used in this library.

\[ \min \sum_{k=1}^l α_k \ s.t. \ c + \sum_{k=1}^l b_{kj} * d_k = v_j \quad \forall j \ -α_k ≤ b_{kj} ≤ α_k \quad \forall k, j \ α_k ≥ 0 \quad \forall k\]

The resulting zonotope has center `c`

and generators `α_k · d_k`

.

Note that the first type of side constraints is vector-based and that the nonnegativity constraints (last type) are not stated explicitly in [1].

[1] *Zonotopes as bounding volumes*. L. J. Guibas, A. T. Nguyen, L. Zhang. SODA 2003.

`LazySets.Approximations._overapproximate_zonotope_cpa`

— Function`_overapproximate_zonotope_cpa(X::LazySet, dir::Type{<:AbstractDirections})`

Overapproximate a polytopic set with a zonotope using Cartesian decomposition.

**Input**

`X`

– polytopic set`dir`

– directions used for the generators

**Output**

A zonotope that overapproximates `X`

.

**Notes**

The algorithm decomposes `X`

into 2D sets and overapproximates those sets with zonotopes, and finally converts the Cartesian product of the sets to a zonotope.

**Algorithm**

The implementation is based on Section 8.2.4 in [1].

[1] Le Guernic, C. *Reachability analysis of hybrid systems with linear continuous dynamics* (Doctoral dissertation). 2009.

## Underapproximations

`LazySets.Approximations.underapproximate`

— Function```
underapproximate(X::S, dirs::AbstractDirections;
[apply_convex_hull]::Bool=false) where {N, S<:LazySet{N}}
```

Compute the underapproximation of a convex set by sampling support vectors.

**Input**

`X`

– convex set`dirs`

– directions`apply_convex_hull`

– (optional, default:`false`

) if`true`

, post-process the support vectors with a convex hull operation

**Output**

The `VPolytope`

obtained by taking the convex hull of support vectors of `X`

along the directions determined by `dirs`

.

**Notes**

Since the support vectors are not always unique, this algorithm may return a strict underapproximation even if the set can be exactly approximated using the given template.

```
underapproximate(X::LazySet, ::Type{<:Hyperrectangle};
solver=nothing) where {N}
```

Underapproximate a convex polygon with a hyperrectangle of maximal area.

**Input**

`X`

– convex polygon`Hyperrectangle`

– type for dispatch`solver`

– (optional; default:`nothing`

) nonlinear solver; if`nothing`

,`default_nln_solver(N)`

will be used

**Output**

A hyperrectangle underapproximation with maximal area.

**Notes**

The implementation only works for 2D sets, but the algorithm can be generalized.

Due to numerical issues, the result may be slightly outside the set.

**Algorithm**

The algorithm is taken from [1, Theorem 17] and solves a convex program (in fact a linear program with nonlinear objective). (The objective is modified to an equivalent version due to solver issues.)

[1] Mehdi Behroozi - *Largest inscribed rectangles in geometric convex sets.* arXiv:1905.13246.

`underapproximate(X::LazySet, ::Type{<:Ball2})`

Compute the largest inscribed Euclidean ball in a set `X`

.

**Input**

`X`

– set`Ball2`

– target type

**Output**

A largest `Ball2`

contained in `X`

.

**Algorithm**

We use `chebyshev_center_radius(X)`

.

## Approximations

`LazySets.Approximations.approximate`

— Function`approximate(R::Rectification; apply_convex_hull::Bool=false)`

Approximate a rectification of a polytopic set with a convex polytope.

**Input**

`R`

– rectification of a polytopic set`apply_convex_hull`

– (optional; default:`false`

) option to remove redundant vertices

**Output**

A polytope in vertex representation (`VPolygon`

in 2D, `VPolytope`

otherwise). There is no guarantee that the result over- or underapproximates `R`

.

**Algorithm**

Let $X$ be the set that is rectified. We compute the vertices of $X$, rectify them, and return the convex hull of the result.

**Notes**

Let $X$ be the set that is rectified and let $p$ and $q$ be two vertices on a facet of $X$. Intuitively, an approximation may occur if the line segment connecting these vertices crosses a coordinate hyperplane and if the line segment connecting the rectified vertices has a different angle.

As a corollary, the approximation is exact for the special cases that the original set is contained in either the positive or negative orthant or is axis-aligned.

## Box Approximations

`LazySets.Approximations.box_approximation`

— Function`box_approximation(S::LazySet{N}) where {N}`

Overapproximate a set by a tight hyperrectangle.

**Input**

`S`

– set

**Output**

A tight hyperrectangle.

**Notes**

An alias for this function is `interval_hull`

.

**Algorithm**

The center and radius of the hyperrectangle are obtained by averaging the low and high coordinates of `S`

computed with the `extrema`

function.

`box_approximation(S::CartesianProductArray{N, <:AbstractHyperrectangle}) where {N}`

Overapproximate the Cartesian product of a finite number of hyperrectangular sets by a tight hyperrectangle.

**Input**

`S`

– Cartesian product of a finite number of hyperrectangular sets

**Output**

A hyperrectangle.

**Algorithm**

This method falls back to the `convert`

method. Since the sets wrapped by the Cartesian product array are hyperrectangles, this can be done without overapproximation.

`box_approximation(S::CartesianProduct{N, <:AbstractHyperrectangle, <:AbstractHyperrectangle}) where {N}`

Overapproximate the Cartesian product of two hyperrectangular sets by a tight hyperrectangle.

**Input**

`S`

– Cartesian product of two hyperrectangular sets

**Output**

A hyperrectangle.

**Algorithm**

This method falls back to the `convert`

method. Since the sets wrapped by the Cartesian product array are hyperrectangles, this can be done without overapproximation.

`box_approximation(lm::LinearMap{N, <:AbstractHyperrectangle}) where {N}`

Return a tight overapproximation of the linear map of a hyperrectangular set using a hyperrectangle.

**Input**

`lm`

– linear map of a hyperrectangular set

**Output**

A hyperrectangle.

**Algorithm**

If `c`

and `r`

denote the center and vector radius of a hyperrectangle `H`

, a tight hyperrectangular overapproximation of `M * H`

is obtained by transforming `c ↦ M*c`

and `r ↦ abs.(M) * r`

, where `abs.(⋅)`

denotes the element-wise absolute-value operator.

`box_approximation(R::Rectification{N}) where {N}`

Overapproximate the rectification of a set by a tight hyperrectangle.

**Input**

`R`

– rectification of a set

**Output**

A hyperrectangle.

**Algorithm**

Box approximation and rectification distribute. We first check whether the wrapped set is empty. If so, we return the empty set. Otherwise, we compute the box approximation of the wrapped set, rectify the resulting box (which is simple), and finally convert the resulting set to a box.

`box_approximation(Z::AbstractZonotope)`

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

**Input**

`Z`

– zonotope

**Output**

A hyperrectangle.

**Algorithm**

This function implements the method in [Section 5.1.2, 1]. A zonotope $Z = ⟨c, G⟩$ can be tightly overapproximated by an axis-aligned 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.

`box_approximation(am::AbstractAffineMap{N, <:AbstractHyperrectangle}) where {N}`

Overapproximate the affine map of a hyperrectangular set by a tight hyperrectangle.

**Input**

`am`

– affine map of a hyperrectangular set

**Output**

A hyperrectangle.

**Algorithm**

If `c`

and `r`

denote the center and vector radius of a hyperrectangle `H`

and `v`

is the translation vector, a tight hyperrectangular overapproximation of `M * H + v`

is obtained by transforming `c ↦ M*c+v`

and `r ↦ abs.(M) * r`

, where `abs.(⋅)`

denotes the element-wise absolute-value operator.

`box_approximation(ch::ConvexHull; [algorithm]::String="box")`

Overapproximate a convex hull with a tight hyperrectangle.

**Input**

`ch`

– convex hull`algorithm`

– (optional; default:`"box"`

) algorithm choice

**Output**

A hyperrectangle.

**Algorithm**

Let `X`

and `Y`

be the two sets of `ch`

. We make use of the following property:

\[\square(CH(X, Y)) = \square\left( X \cup Y \right) = \square\left( \square(X) \cup \square(Y) \right)\]

If `algorithm == "extrema"`

, we compute the low and high coordinates of `X`

and `Y`

via `extrema`

.

If `algorithm == "box"`

, we instead compute the box approximations of `X`

and `Y`

via `box_approximation`

.

In both cases we then take the box approximation of the result.

The `"extrema"`

algorithm is more efficient if `extrema`

is efficient because it does not need to allocate the intermediate hyperrectangles.

`box_approximation(ms::MinkowskiSum)`

Overapproximate the Minkowski sum of two sets with a tight hyperrectangle.

**Input**

`ms`

– Minkowski sum

**Output**

A hyperrectangle.

**Algorithm**

The box approximation distributes over the Minkowski sum:

\[\square(X \oplus Y) = \square(X) \oplus \square(Y)\]

It suffices to compute the box approximation of each summand and then take the concrete Minkowski sum for hyperrectangles.

`LazySets.Approximations.interval_hull`

— Function`interval_hull`

Alias for `box_approximation`

.

`LazySets.Approximations.symmetric_interval_hull`

— Function`symmetric_interval_hull(S::LazySet{N}) where {N}`

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

**Input**

`S`

– set

**Output**

A tight hyperrectangle that is centrally symmetric wrt. the origin.

**Notes**

An alias for this function is `box_approximation_symmetric`

.

**Algorithm**

The center of the box is the origin, and the radius is obtained via the `extrema`

function.

`LazySets.Approximations.box_approximation_symmetric`

— Function`box_approximation_symmetric`

Alias for `symmetric_interval_hull`

.

`LazySets.Approximations.ballinf_approximation`

— Function`ballinf_approximation(S::LazySet)`

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

**Input**

`S`

– set

**Output**

A tight ball in the infinity norm.

**Algorithm**

The center and radius of the ball are obtained by averaging the low and high coordinates of `S`

computed with the `extrema`

function.

## Iterative refinement

`LazySets.Approximations.overapproximate_hausdorff`

— Function`overapproximate_hausdorff(X::S, ε::Real) where {N<:AbstractFloat, S<:LazySet{N}}`

Return an ε-close overapproximation of the given 2D convex set (in terms of the Hausdorff distance) in the form of a polygon in constraint representation.

**Input**

`X`

– 2D convex set`ε`

– error bound

**Output**

A polygon in constraint representation.

`LazySets.Approximations.LocalApproximation`

— Type`LocalApproximation{N, VN<:AbstractVector{N}}`

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`

– flag stating whether this approximation is refinable`err`

– error upper bound

**Notes**

The criteria for being refinable are determined in `new_approx`

.

`LazySets.Approximations.PolygonalOverapproximation`

— Type`PolygonalOverapproximation{N, SN<:LazySet{N}, VN<:AbstractVector{N}}`

Type that represents a polygonal overapproximation of a convex set.

**Fields**

`S`

– convex set`approx_stack`

– stack of local approximations that still need to be examined`constraints`

– vector of half-spaces that are already finalized (i.e., they satisfy the given error bound)

`LazySets.Approximations.new_approx`

— Function```
new_approx(S::LazySet, p1::VN, d1::VN,
p2::VN, d2::VN) where {N<:AbstractFloat, VN<:AbstractVector{N}}
```

Create a `LocalApproximation`

instance for the given excerpt of a polygonal overapproximation.

**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.

`LazySets.Approximations.addapproximation!`

— Function```
addapproximation!(Ω::PolygonalOverapproximation, p1::VN, d1::VN,
p2::VN, d2::VN) where {N, VN<:AbstractVector{N}}
```

**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.

`LazySets.Approximations.refine`

— Method`refine(approx::LocalApproximation, S::LazySet)`

Refine a given local approximation of the polygonal overapproximation 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.

`LazySets.Approximations.tohrep`

— Method`tohrep(Ω::PolygonalOverapproximation)`

Convert a polygonal overapproximation into a polygon in constraint representation.

**Input**

`Ω`

– polygonal overapproximation of a convex set

**Output**

A polygon in constraint representation.

**Algorithm**

Internally, the constraints of `Ω`

are already sorted.

`Base.convert`

— Method`convert(::Type{HalfSpace}, approx::LocalApproximation)`

Convert a local approximation to a half-space.

**Input**

`approx`

– local approximation

**Output**

A half-space.

## Template directions

`LazySets.Approximations.AbstractDirections`

— Type`AbstractDirections{N, VN}`

Abstract type for representations of direction vectors.

**Notes**

This type is parameterized by `N`

and `VN`

, where:

`N`

stands for the numeric type`VN`

stands for the vector type with coefficients of type`N`

Each implementing subtype is an iterator over a set of directions. For that they implement the standard iterator methods from `Base`

, namely `Base.length`

(returns the number of directions) and `Base.iterate`

. Moreover, the following methods should be implemented:

`dim`

– return the ambient dimension of the vectors`eltype`

– return the type of each vector

Optionally, subtypes may implement:

`isbounding`

– (defaults to`false`

) return`true`

if an overapproximation with the direction vectors results in a bounded set, given a bounded input set, and`false`

otherwise`isnormalized`

– (defaults to`false`

) is`true`

if each direction vector has norm one w.r.t. the usual vector 2-norm

`LazySets.Approximations.isbounding`

— Function```
isbounding(ad::AbstractDirections)
isbounding(ad::Type{<:AbstractDirections})
```

Check whether an overapproximation with a set of direction vectors results in a bounded set, given a bounded input set.

**Input**

`ad`

– direction vectors or a subtype of`AbstractDirections`

**Output**

Given a bounded set $X$, we can construct an outer polyhedral approximation of $X$ by using the direction vectors `ad`

as normal vectors of the facets. If this function returns `true`

, then the result is again guaranteed to be a bounded set (i.e., a polytope). Note that the result does not depend on the specific shape of $X$, as long as $X$ is bounded.

**Notes**

By default, this function returns `false`

in order to be conservative. Custom subtypes of `AbstractDirections`

should hence add a method for this function.

The function can be applied to an instance of an `AbstractDirections`

subtype or to the subtype itself. By default, the check on the instance falls back to the check on the subtype.

`LazySets.Approximations.isnormalized`

— Function```
isnormalized(ad::AbstractDirections)
isnormalized(ad::Type{<:AbstractDirections})
```

Check whether the given direction vectors are normalized with respect to the 2-norm.

**Input**

`ad`

– direction vectors or a subtype of`AbstractDirections`

**Output**

`true`

if the 2-norm of each element in `ad`

is one and `false`

otherwise.

**Notes**

By default, this function returns `false`

in order to be conservative. Custom subtypes of `AbstractDirections`

should hence add a method for this function.

The function can be applied to an instance of an `AbstractDirections`

subtype or to the subtype itself. By default, the check on the instance falls back to the check on the subtype.

`LazySets.project`

— Method```
project(S::LazySet,
block::AbstractVector{Int},
directions::Type{<:AbstractDirections},
[n]::Int;
[kwargs...]
)
```

Project a high-dimensional set to a given block using direction vectors.

**Input**

`S`

– set`block`

– block structure - a vector with the dimensions of interest`directions`

– direction vectors`n`

– (optional, default:`dim(S)`

) ambient dimension of the set`S`

**Output**

The polyhedral overapproximation of the projection of `S`

in the given directions.

`LazySets.Approximations.BoxDirections`

— Type`BoxDirections{N, VN} <: AbstractDirections{N, VN}`

Box directions representation.

**Fields**

`n`

– dimension

**Notes**

Box directions can be seen as the vectors where only one entry is ±1, and all other entries are 0. In dimension $n$, there are $2n$ such directions.

The default vector representation used in this template is a `ReachabilityBase.Arrays.SingleEntryVector`

, although other implementations can be used such as a regular `Vector`

and a `SparseVector`

.

**Examples**

The template can be constructed by passing the dimension. For example, in dimension two:

```
julia> dirs = BoxDirections(2)
BoxDirections{Float64, ReachabilityBase.Arrays.SingleEntryVector{Float64}}(2)
julia> length(dirs)
4
```

By default, each direction is represented as a `SingleEntryVector`

, i.e., a vector with only one non-zero element,

```
julia> eltype(dirs)
ReachabilityBase.Arrays.SingleEntryVector{Float64}
```

In two dimensions, the directions defined by `BoxDirections`

are normal to the facets of a box.

```
julia> collect(dirs)
4-element Vector{ReachabilityBase.Arrays.SingleEntryVector{Float64}}:
[1.0, 0.0]
[0.0, 1.0]
[0.0, -1.0]
[-1.0, 0.0]
```

The numeric type can be specified as well:

```
julia> BoxDirections{Rational{Int}}(10)
BoxDirections{Rational{Int64}, ReachabilityBase.Arrays.SingleEntryVector{Rational{Int64}}}(10)
julia> length(ans)
20
```

`LazySets.Approximations.DiagDirections`

— Type`DiagDirections{N, VN} <: AbstractDirections{N, VN}`

Diagonal directions representation.

**Fields**

`n`

– dimension

**Notes**

Diagonal directions are vectors where all entries are ±1. In dimension $n$, there are in total $2^n$ such directions.

**Examples**

The template can be constructed by passing the dimension. For example, in dimension two:

```
julia> dirs = DiagDirections(2)
DiagDirections{Float64, Vector{Float64}}(2)
julia> length(dirs) # number of directions
4
```

By default, each direction is represented as a regular `Vector`

:

```
julia> eltype(dirs)
Vector{Float64} (alias for Array{Float64, 1})
```

In two dimensions, the directions defined by `DiagDirections`

are normal to the facets of a ball in the 1-norm.

```
julia> collect(dirs)
4-element Vector{Vector{Float64}}:
[1.0, 1.0]
[-1.0, 1.0]
[1.0, -1.0]
[-1.0, -1.0]
```

The numeric type can be specified as well:

```
julia> DiagDirections{Rational{Int}}(10)
DiagDirections{Rational{Int64}, Vector{Rational{Int64}}}(10)
julia> length(ans)
1024
```

`LazySets.Approximations.OctDirections`

— Type`OctDirections{N, VN} <: AbstractDirections{N, VN}`

Octagon directions 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$. In dimension $n$, there are $2n^2$ such directions.

**Examples**

The template can be constructed by passing the dimension. For example, in dimension two:

```
julia> dirs = OctDirections(2)
OctDirections{Float64, SparseArrays.SparseVector{Float64, Int64}}(2)
julia> length(dirs) # number of directions
8
```

By default, the directions are represented as sparse vectors:

```
julia> eltype(dirs)
SparseArrays.SparseVector{Float64, Int64}
```

In two dimensions, the directions are normal to the facets of an octagon.

```
julia> first(dirs)
2-element SparseArrays.SparseVector{Float64, Int64} with 2 stored entries:
[1] = 1.0
[2] = 1.0
julia> Vector.(collect(dirs))
8-element Vector{Vector{Float64}}:
[1.0, 1.0]
[1.0, -1.0]
[-1.0, 1.0]
[-1.0, -1.0]
[1.0, 0.0]
[0.0, 1.0]
[0.0, -1.0]
[-1.0, 0.0]
```

The numeric type can be specified as well:

```
julia> OctDirections{Rational{Int}}(10)
OctDirections{Rational{Int64}, SparseArrays.SparseVector{Rational{Int64}, Int64}}(10)
julia> length(ans)
200
```

`LazySets.Approximations.BoxDiagDirections`

— Type`BoxDiagDirections{N, VN} <: AbstractDirections{N, VN}`

Box-diagonal directions 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. In dimension $n$, there are in total $2^n + 2n$ such directions.

**Examples**

The template can be constructed by passing the dimension. For example, in dimension two:

```
julia> dirs = BoxDiagDirections(2)
BoxDiagDirections{Float64, Vector{Float64}}(2)
julia> length(dirs) # number of directions
8
```

By default, each direction is represented as a regular vector:

```
julia> eltype(dirs)
Vector{Float64} (alias for Array{Float64, 1})
```

In two dimensions, the directions are normal to the facets of an octagon, i.e., the template coincides with `OctDirections`

.

```
julia> collect(dirs)
8-element Vector{Vector{Float64}}:
[1.0, 1.0]
[-1.0, 1.0]
[1.0, -1.0]
[-1.0, -1.0]
[1.0, 0.0]
[0.0, 1.0]
[0.0, -1.0]
[-1.0, 0.0]
```

The numeric type can be specified as well:

```
julia> BoxDiagDirections{Rational{Int}}(10)
BoxDiagDirections{Rational{Int64}, Vector{Rational{Int64}}}(10)
julia> length(ans)
1044
```

`LazySets.Approximations.PolarDirections`

— Type`PolarDirections{N<:AbstractFloat, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}`

Polar directions representation.

**Fields**

`Nφ`

– length of the partition of the polar angle`stack`

– list of computed directions

**Notes**

The `PolarDirections`

constructor computes a sample of the unit sphere in $\mathbb{R}^2$, which is parameterized by the polar angle $φ ∈ Dφ := [0, 2π]$; see the Wikipedia entry on the polar coordinate system for details. The resulting directions are stored in `stack`

.

The integer argument $Nφ$ defines how many samples of $Dφ$ are taken. The Cartesian components of each direction are obtained with

\[[cos(φᵢ), sin(φᵢ)].\]

**Examples**

The integer passed as an argument is used to discretize $φ$:

```
julia> pd = PolarDirections(2);
julia> pd.stack
2-element Vector{Vector{Float64}}:
[1.0, 0.0]
[-1.0, 1.2246467991473532e-16]
julia> length(pd)
2
```

`LazySets.Approximations.SphericalDirections`

— Type`SphericalDirections{N<:AbstractFloat, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}`

Spherical directions representation.

**Fields**

`Nθ`

– length of the partition of the azimuthal angle`Nφ`

– length of the partition of the polar angle`stack`

– list of computed directions

**Notes**

The `SphericalDirections`

constructor provides a sample of the unit sphere in $\mathbb{R}^3$, which is parameterized by the azimuthal and polar angles $θ ∈ Dθ := [0, π]$ and $φ ∈ Dφ := [0, 2π]$ respectively; see the Wikipedia entry on the spherical coordinate system for details.

The integer arguments $Nθ$ and $Nφ$ define how many samples along the domains $Dθ$ and $Dφ$ are respectively taken. The Cartesian components of each direction are obtained with

\[[sin(θᵢ)*cos(φᵢ), sin(θᵢ)*sin(φᵢ), cos(θᵢ)].\]

The north and south poles are treated separately so that those points are not considered more than once.

**Examples**

The template can be built in different ways. If you pass only one integer, the same value is used to discretize both $θ$ and $φ$:

```
julia> sd = SphericalDirections(3);
julia> sd.Nθ, sd.Nφ
(3, 3)
julia> length(sd)
4
```

Pass two integers to control the discretization in $θ$ and in $φ$ separately:

```
julia> sd = SphericalDirections(4, 5);
julia> length(sd)
10
julia> sd = SphericalDirections(4, 8);
julia> length(sd)
16
```

`LazySets.Approximations.CustomDirections`

— Type`CustomDirections{N, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}`

User-defined direction vectors.

**Fields**

`directions`

– list of direction vectors`n`

– (optional; default: computed from`directions`

) dimension`check_boundedness`

– (optional; default:`true`

) flag to check boundedness`check_normalization`

– (optional; default:`true`

) flag to check whether all directions are normalized

**Notes**

This struct is a wrapper for a list of user-defined directions. There are fields for the list of directions, their dimension, and (boolean) cache fields for the boundedness and normalization properties. The latter are checked by default upon construction.

To check boundedness, we construct the polyhedron with constraints $d·x <= 1$ for each direction $d$ and check if this set is bounded. (Note that the bound $1$ is arbitrary and that this set may be empty, which however implies boundedness.)

The dimension will also be determined automatically, unless the empty vector is passed (in which case the optional argument `n`

needs to be specified).

**Examples**

Create a template with box directions in dimension two:

```
julia> dirs = CustomDirections([[1.0, 0.0], [-1.0, 0.0], [0.0, 1.0], [0.0, -1.0]]);
julia> dirs.directions
4-element Vector{Vector{Float64}}:
[1.0, 0.0]
[-1.0, 0.0]
[0.0, 1.0]
[0.0, -1.0]
julia> LazySets.Approximations.isbounding(dirs)
true
julia> LazySets.Approximations.isnormalized(dirs)
true
```

See also `overapproximate(X::LazySet, dir::AbstractDirections)`

.

## Hausdorff distance

`LazySets.Approximations.hausdorff_distance`

— Function```
hausdorff_distance(X::LazySet{N}, Y::LazySet{N}; [p]::N=N(Inf),
[ε]=N(1e-3)) where {N}
```

Compute the Hausdorff distance between two convex sets up to a given threshold.

**Input**

`X`

– convex set`Y`

– convex set`p`

– (optional, default:`Inf`

) norm parameter of the Hausdorff distance`ε`

– (optional, default:`1e-3`

) precision threshold; the true Hausdorff distance is allowed to diverge from the result by at most this value

**Output**

A value from the $ε$-neighborhood of the Hausdorff distance between $X$ and $Y$.

**Notes**

Given a $p$-norm, the Hausdorff distance $d_H^p(X, Y)$ between sets $X$ and $Y$ is defined as follows:

\[ d_H^p(X, Y) = \inf\{δ ≥ 0 \mid Y ⊆ X ⊕ δ 𝐵_p^n \text{ and } X ⊆ Y ⊕ δ 𝐵_p^n\}\]

Here $𝐵_p^n$ is the $n$-dimensional unit ball in the $p$-norm.

The implementation may internally rely on the support function of $X$ and $Y$; hence any imprecision in the implementation of the support function may affect the result. At the time of writing, the only convex set type with imprecise support function is the lazy `Intersection`

.

**Algorithm**

We perform binary search for bounding the Hausdorff distance in an interval $[l, u]$, where initially $l$ is $0$ and $u$ is described below. The binary search terminates when $u - l ≤ ε$, i.e., the interval becomes sufficiently small.

To find an upper bound $u$, we start with the heuristics of taking the biggest distance in the axis-parallel directions. As long as this bound does not work, we increase the bound by $2$.

Given a value $δ$, to check whether the sets are within Hausdorff distance $δ$, we simply check the inclusions given above, where on the right-hand side we use a lazy `Bloating`

.