# Polyhedra (AbstractPolyhedron)

A polyhedron has finitely many facets (*H-representation*) and is not necessarily bounded.

`LazySets.AbstractPolyhedron`

— Type`AbstractPolyhedron{N} <: ConvexSet{N}`

Abstract type for compact convex polyhedral sets.

**Notes**

Every concrete `AbstractPolyhedron`

must define the following functions:

`constraints_list(::AbstractPolyhedron{N})`

– return a list of all facet constraints

Polyhedra are defined as the intersection of a finite number of closed half-spaces. As such, polyhedra are closed and convex but not necessarily bounded. Bounded polyhedra are called *polytopes* (see `AbstractPolytope`

).

The subtypes of `AbstractPolyhedron`

(including abstract interfaces):

```
julia> subtypes(AbstractPolyhedron)
8-element Vector{Any}:
AbstractPolytope
HPolyhedron
HalfSpace
Hyperplane
Line
Line2D
Star
Universe
```

This interface defines the following functions:

`Base.:∈`

— Method`∈(x::AbstractVector, P::AbstractPolyhedron)`

Check whether a given point is contained in a polyhedron.

**Input**

`x`

– point/vector`P`

– polyhedron

**Output**

`true`

iff $x ∈ P$.

**Algorithm**

This implementation checks if the point lies inside each defining half-space.

`LazySets.isuniversal`

— Method`isuniversal(P::AbstractPolyhedron{N}, [witness]::Bool=false) where {N}`

Check whether a polyhedron is universal.

**Input**

`P`

– polyhedron`witness`

– (optional, default:`false`

) compute a witness if activated

**Output**

- If
`witness`

option is deactivated:`true`

iff $P$ is universal - If
`witness`

option is activated:`(true, [])`

iff $P$ is universal`(false, v)`

iff $P$ is not universal and $v ∉ P$

**Algorithm**

`P`

is universal iff it has no constraints.

A witness is produced using `isuniversal(H)`

where `H`

is the first linear constraint of `P`

.

`LazySets.constrained_dimensions`

— Method`constrained_dimensions(P::AbstractPolyhedron)`

Return the indices in which a polyhedron is constrained.

**Input**

`P`

– polyhedron

**Output**

A vector of ascending indices `i`

such that the polyhedron is constrained in dimension `i`

.

**Examples**

A 2D polyhedron with constraint $x1 ≥ 0$ is constrained in dimension 1 only.

`LazySets.an_element`

— Method```
an_element(P::AbstractPolyhedron{N};
[solver]=default_lp_solver(N)) where {N}
```

Return some element of a polyhedron.

**Input**

`P`

– polyhedron`solver`

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

) LP solver

**Output**

An element of the polyhedron, or an error if the polyhedron is empty.

**Algorithm**

An element is obtained by solving a feasibility linear program.

`an_element(U::Universe{N}) where {N}`

Return some element of a universe.

**Input**

`U`

– universe

**Output**

The origin.

`LazySets.isbounded`

— Method`isbounded(P::AbstractPolyhedron{N}; [solver]=default_lp_solver(N)) where {N}`

Check whether a polyhedron is bounded.

**Input**

`P`

– polyhedron`solver`

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

) the backend used to solve the linear program

**Output**

`true`

iff the polyhedron is bounded

**Algorithm**

We first check if the polyhedron has more than `dim(P)`

constraints, which is a necessary condition for boundedness.

If so, we check boundedness via `_isbounded_stiemke`

.

`LazySets.vertices_list`

— Method`vertices_list(P::AbstractPolyhedron; check_boundedness::Bool=true)`

Return the list of vertices of a polyhedron in constraint representation.

**Input**

`P`

– polyhedron in constraint representation`check_boundedness`

– (optional, default:`true`

) if`true`

, check whether the polyhedron is bounded

**Output**

The list of vertices of `P`

, or an error if `P`

is unbounded.

**Notes**

This function throws an error if the polyhedron is unbounded. Otherwise, the polyhedron is converted to an `HPolytope`

and its list of vertices is computed.

**Examples**

```
julia> P = HPolyhedron([HalfSpace([1.0, 0.0], 1.0),
HalfSpace([0.0, 1.0], 1.0),
HalfSpace([-1.0, 0.0], 1.0),
HalfSpace([0.0, -1.0], 1.0)]);
julia> length(vertices_list(P))
4
```

`LazySets.project`

— Method```
project(P::AbstractPolyhedron{N}, block::AbstractVector{Int};
[kwargs...]) where {N}
```

Concrete projection of a polyhedral set.

**Input**

`P`

– set`block`

– block structure, a vector with the dimensions of interest

**Output**

A polyhedron representing the projection of `P`

on the dimensions specified by `block`

. If `P`

was bounded, the result is an `HPolytope`

; otherwise the result is an `HPolyhedron`

. Note that there are more specific methods for specific input types, which give a different output type; e.g., projecting a `Ball1`

results in a `Ball1`

.

**Algorithm**

- We first try to exploit the special case where each of the constraints of
`P`

and`block`

are*compatible*, which is one of the two cases described below. Let`c`

be a constraint of`P`

and let $D_c$ and $D_b$ be the set of dimensions in which`c`

resp.`block`

are constrained.- If $D_c ⊆ D_b$, then one can project the normal vector of
`c`

. - If $D_c ∩ D_b = ∅$, then the constraint becomes redundant.

- If $D_c ⊆ D_b$, then one can project the normal vector of
- In the general case, we compute the concrete linear map of the projection matrix associated to the given block structure.

**Examples**

Consider the four-dimensional cross-polytope (unit ball in the 1-norm):

`julia> P = convert(HPolytope, Ball1(zeros(4), 1.0));`

All dimensions are constrained, and computing the (trivial) projection on the whole space behaves as expected:

```
julia> constrained_dimensions(P)
4-element Vector{Int64}:
1
2
3
4
julia> project(P, [1, 2, 3, 4]) == P
true
```

Each constraint of the cross polytope is constrained in all dimensions.

Now let us take a ball in the infinity norm and remove some constraints:

```
julia> B = BallInf(zeros(4), 1.0);
julia> c = constraints_list(B)[1:2]
2-element Vector{HalfSpace{Float64, ReachabilityBase.Arrays.SingleEntryVector{Float64}}}:
HalfSpace{Float64, ReachabilityBase.Arrays.SingleEntryVector{Float64}}([1.0, 0.0, 0.0, 0.0], 1.0)
HalfSpace{Float64, ReachabilityBase.Arrays.SingleEntryVector{Float64}}([0.0, 1.0, 0.0, 0.0], 1.0)
julia> P = HPolyhedron(c);
julia> constrained_dimensions(P)
2-element Vector{Int64}:
1
2
```

Finally, we take the concrete projection onto variables `1`

and `2`

:

```
julia> project(P, [1, 2]) |> constraints_list
2-element Vector{HalfSpace{Float64, Vector{Float64}}}:
HalfSpace{Float64, Vector{Float64}}([1.0, 0.0], 1.0)
HalfSpace{Float64, Vector{Float64}}([0.0, 1.0], 1.0)
```

`LazySets._linear_map_polyhedron`

— Function```
_linear_map_polyhedron(M::AbstractMatrix{NM},
P::LazySet{NP};
[algorithm]::Union{String, Nothing}=nothing,
[check_invertibility]::Bool=true,
[cond_tol]::Number=DEFAULT_COND_TOL,
[inverse]::Union{AbstractMatrix{N}, Nothing}=nothing,
[backend]=nothing,
[elimination_method]=nothing) where {NM, NP}
```

Concrete linear map of a polyhedral set.

**Input**

`M`

– matrix`P`

– polyhedral set`algorithm`

– (optional; default:`nothing`

) algorithm to be used; for the description see the Algorithm section below; possible choices are:`"inverse"`

, alias:`"inv"`

`"inverse_right"`

, alias:`"inv_right"`

`"elimination"`

, alias:`"elim"`

`"lift"`

`"vrep"`

`"vrep_chull"`

`check_invertibility`

– (optional, default:`true`

) if`true`

check whether the given matrix`M`

is invertible; set to`false`

only if you know that`M`

is invertible`cond_tol`

– (optional; default:`DEFAULT_COND_TOL`

) tolerance of matrix condition (used to check whether the matrix is invertible)`inverse`

– (optional; default:`nothing`

) matrix inverse`M⁻¹`

; use this option if you have already computed the inverse matrix of`M`

`backend`

– (optional: default:`nothing`

) polyhedra backend`elimination_method`

– (optional: default:`nothing`

) elimination method for the`"elimination"`

algorithm

**Output**

The type of the result is "as close as possible" to the the type of `P`

. Let `(m, n)`

be the size of `M`

, where `m ≠ n`

is allowed for rectangular maps.

To fix the type of the output to something different than the default value, consider post-processing the result of this function with a call to a suitable `convert`

method.

In particular, the output depends on the type of `P`

, on `m`

, and the algorithm that was used:

If the vertex-based approach was used:

- If
`P`

is a`VPolygon`

and`m = 2`

then the output is a`VPolygon`

. - If
`P`

is a`VPolytope`

then the output is a`VPolytope`

. - Otherwise the output is an
`Interval`

if`m = 1`

, a`VPolygon`

if`m = 2`

, and a`VPolytope`

in other cases.

- If
If the invertibility criterion was used:

- The types of
`HalfSpace`

,`Hyperplane`

,`Line2D`

, and subtypes of`AbstractHPolygon`

are preserved. - If
`P`

is an`AbstractPolytope`

, then the output is an`Interval`

if`m = 1`

, an`HPolygon`

if`m = 2`

, and an`HPolytope`

in other cases. - Otherwise the output is an
`HPolyhedron`

.

- The types of

**Notes**

Since the different linear-map algorithms work at the level of constraints, this method uses dispatch on two stages: once the algorithm has been defined, first the helper methods `_linear_map_hrep_helper`

(resp. `_linear_map_vrep`

) are invoked, which dispatch on the set type. Then, each helper method calls the concrete implementation of `_linear_map_hrep`

, which dispatches on the algorithm, and returns a list of constraints.

To simplify working with different algorithms and options, the types `<: AbstractLinearMapAlgorithm`

are used. These types are singleton type or types that carry only the key data for the given algorithm, such as the matrix inverse or the polyhedra backend.

New subtypes of the `AbstractPolyhedron`

interface may define their own helper methods `_linear_map_vrep`

(respectively `_linear_map_hrep_helper`

) for special handling of the constraints returned by the implementations of `_linear_map_hrep`

; otherwise the fallback implementation for `AbstractPolyhedron`

is used, which instantiates an `HPolyhedron`

.

**Algorithm**

This method mainly implements several approaches for the linear map: inverse, right inverse, transformation to vertex representation, variable elimination, and variable lifting. Depending on the properties of `M`

and `P`

, one algorithm may be preferable over the other. Details on the algorithms are given in the following subsections.

Otherwise, if the algorithm argument is not specified, a default option is chosen based on heuristics on the types and values of `M`

and `P`

:

- If the
`"inverse"`

algorithm applies, it is used. - Otherwise, if the
`"inverse_right"`

algorithm applies, it is used. - Otherwise, if the
`"lift"`

algorithm applies, it is used. - Otherwise, the
`"elimination"`

algorithm is used.

Note that the algorithms `"inverse"`

and `"inverse_right"`

do not require the external library `Polyhedra`

. However, the fallback method `"elimination"`

requires `Polyhedra`

as well as the library `CDDLib`

.

The optional keyword arguments `inverse`

and `check_invertibility`

modify the default behavior:

- If an inverse matrix is passed in
`inverse`

, the given algorithm is applied, and if none is given, either`"inverse"`

or`"inverse_right"`

is applied (in that order of preference). - If
`check_invertibility`

is set to`false`

, the given algorithm is applied, and if none is given, either`"inverse"`

or`"inverse_right"`

is applied (in that order of preference).

**Inverse**

This algorithm is invoked with the keyword argument `algorithm="inverse"`

(or `algorithm="inv"`

). The algorithm requires that `M`

is invertible, square, and dense. If you know a priori that `M`

is invertible, set the flag `check_invertibility=false`

, such that no extra checks are done. Otherwise, we check the sufficient condition that the condition number of `M`

is not too high. The threshold for the condition number can be modified from its default value, `DEFAULT_COND_TOL`

, by passing a custom `cond_tol`

.

The algorithm is described next. Assuming that the matrix $M$ is invertible (which we check via a sufficient condition,), $y = M x$ implies $x = \text{inv}(M) y$ and we can transform the polyhedron $A x ≤ b$ to the polyhedron $A \text{inv}(M) y ≤ b$.

If the dense condition on `M`

is not satisfied, there are two suggested workarounds: either transform to a dense matrix, i.e., calling `linear_map`

with `Matrix(M)`

, or use the `"inverse_right"`

algorithm, which does not compute the inverse matrix explicitly, but uses a polyalgorithm; see the documentation of `?`

for details.

**Inverse-right**

This algorithm is invoked with the keyword argument `algorithm="inverse_right"`

(or `algorithm="inv_right"`

). This algorithm applies to square and invertible matrices `M`

. The idea is essentially the same as for the `"inverse"`

algorithm; the difference is that in `"inverse"`

the full matrix inverse is computed, and in `"inverse_right"`

only the left division on the normal vectors is used. In particular, `"inverse_right"`

is good as a workaround when `M`

is sparse (since the `inv`

function is not available for sparse matrices).

**Elimination**

This algorithm is invoked with the keyword argument `algorithm = "elimination"`

(or `algorithm = "elim"`

). The algorithm applies to any matrix `M`

(invertible or not), and any polyhedron `P`

(bounded or not).

The idea is described next. If `P : Ax <= b`

and `y = Mx`

denote the polyhedron and the linear map, respectively, we consider the vector `z = [y, x]`

, write the given equalities and the inequalities, and then eliminate the last x variables (there are `length(x)`

in total) using a call to `Polyhedra.eliminate`

to a backend library that can do variable elimination (typically `CDDLib`

with the `BlockElimination()`

algorithm). In this way we have eliminated the "old" variables `x`

and kept the "new" or transformed variables "y".

The default elimination method is block elimination. For possible options we refer to the documentation of Polyhedra, projection/elimination.

**Lift**

This algorithm is invoked with the keyword argument `algorithm="lift"`

. The algorithm applies if `M`

is rectangular of size `m × n`

with `m > n`

and full rank (i.e., of rank `n`

).

The idea is to embed the polyhedron into the `m`

-dimensional space by appending zeros, i.e. extending all constraints of `P`

to `m`

dimensions, and constraining the last `m - n`

dimensions to `0`

. The resulting matrix is extended to an invertible `m × m`

matrix, and the algorithm using the inverse of the linear map is applied. For technical details of extending `M`

to a higher-dimensional invertible matrix, see `ReachabilityBase.Arrays.extend`

.

**Vertex representation**

This algorithm is invoked with the keyword argument `algorithm="vrep"`

(or `algorithm="vrep_chull"`

). If the polyhedron is bounded, the idea is to convert it to its vertex representation and apply the linear map to each vertex.

The returned set is a polytope in vertex representation. Note that conversion of the result back to half-space representation is not computed by default, since this may be costly. If you use this algorithm and still want to convert back to half-space representation, apply `tohrep`

to the result of this method.

`LazySets._isbounded_stiemke`

— Function```
_isbounded_stiemke(constraints::AbstractVector{<:HalfSpace{N}};
solver=LazySets.default_lp_solver(N),
check_nonempty::Bool=true) where {N}
```

Check whether a list of constraints is bounded using Stiemke's theorem of alternatives.

**Input**

`constraints`

– list of constraints`backend`

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

) the backend used to solve the linear program`check_nonempty`

– (optional, default:`true`

) if`true`

, check the precondition to this algorithm that`P`

is non-empty

**Output**

`true`

iff the list of constraints is bounded.

**Notes**

The list of constraints represents a polyhedron.

The algorithm calls `isempty`

to check whether the polyhedron is empty. This computation can be avoided using the `check_nonempty`

flag.

**Algorithm**

The algorithm is based on Stiemke's theorem of alternatives, see, e.g., [1].

Let the polyhedron $P$ be given in constraint form $Ax ≤ b$. We assume that the polyhedron is non-empty.

Proposition 1. If $\ker(A)≠\{0\}$, then $P$ is unbounded.

Proposition 2. Assume that $ker(A)={0}$ and $P$ is non-empty. Then $P$ is bounded if and only if the following linear program admits a feasible solution: $\min∥y∥_1$ subject to $A^Ty=0$ and $y≥1$.

[1] Mangasarian, Olvi L. *Nonlinear programming.* Society for Industrial and Applied Mathematics, 1994.

Some common functions to work with linear constraints:

`LazySets.constraints_list`

— Method`constraints_list(A::AbstractMatrix{N}, b::AbstractVector)`

Convert a matrix-vector representation to a linear-constraint representation.

**Input**

`A`

– matrix`b`

– vector

**Output**

A list of linear constraints.

`LazySets.tosimplehrep`

— Method```
tosimplehrep(constraints::AbstractVector{LC};
[n]::Int=0) where {N, LC<:HalfSpace{N}}
```

Return the simple H-representation $Ax ≤ b$ from a list of linear constraints.

**Input**

`constraints`

– a list of linear constraints`n`

– (optional; default:`0`

) dimension of the constraints

**Output**

The tuple `(A, b)`

where `A`

is the matrix of normal directions and `b`

is the vector of offsets.

**Notes**

The parameter `n`

can be used to create a matrix with no constraints but a non-zero dimension.

`LazySets.remove_redundant_constraints`

— Method```
remove_redundant_constraints(constraints::AbstractVector{S};
backend=nothing) where {S<:HalfSpace}
```

Remove the redundant constraints of a given list of linear constraints.

**Input**

`constraints`

– list of constraints`backend`

– (optional, default:`nothing`

) the backend used to solve the linear program

**Output**

The list of constraints with the redundant ones removed, or an empty set if the constraints are infeasible.

**Notes**

If `backend`

is `nothing`

, it defaults to `default_lp_solver(N)`

.

**Algorithm**

See `remove_redundant_constraints!(::AbstractVector{<:HalfSpace})`

for details.

`LazySets.remove_redundant_constraints!`

— Method```
remove_redundant_constraints!(constraints::AbstractVector{S};
[backend]=nothing) where {S<:HalfSpace}
```

Remove the redundant constraints of a given list of linear constraints; the list is updated in-place.

**Input**

`constraints`

– list of constraints`backend`

– (optional, default:`nothing`

) the backend used to solve the linear program

**Output**

`true`

if the removal was successful and the list of constraints `constraints`

is modified by removing the redundant constraints, and `false`

only if the constraints are infeasible.

**Notes**

Note that the result may be `true`

even if the constraints are infeasible. For example, $x ≤ 0 && x ≥ 1$ will return `true`

without removing any constraint. To check if the constraints are infeasible, use `isempty(HPolyhedron(constraints))`

.

If `backend`

is `nothing`

, it defaults to `default_lp_solver(N)`

.

**Algorithm**

If there are `m`

constraints in `n`

dimensions, this function checks one by one if each of the `m`

constraints is implied by the remaining ones.

To check if the `k`

-th constraint is redundant, an LP is formulated using the constraints that have not yet been removed. If, at an intermediate step, it is detected that a subgroup of the constraints is infeasible, this function returns `false`

. If the calculation finished successfully, this function returns `true`

.

For details, see Fukuda's Polyhedra FAQ.

Plotting (bounded) polyhedra is available, too:

`LazySets.plot_recipe`

— Method`plot_recipe(P::AbstractPolyhedron{N}, [ε]=zero(N)) where {N}`

Convert a (bounded) polyhedron to a pair `(x, y)`

of points for plotting.

**Input**

`P`

– bounded polyhedron`ε`

– (optional, default:`0`

) ignored, used for dispatch

**Output**

A pair `(x, y)`

of points that can be plotted, where `x`

is the vector of x-coordinates and `y`

is the vector of y-coordinates.

**Algorithm**

We first assert that `P`

is bounded (i.e., that `P`

is a polytope).

One-dimensional polytopes are converted to an `Interval`

. Three-dimensional or higher-dimensional polytopes are not supported.

For two-dimensional polytopes (i.e., polygons) we compute their set of vertices using `vertices_list`

and then plot the convex hull of these vertices.