Polyhedra (AbstractPolyhedron)

AbstractPolyhedron{N} <: ConvexSet{N}

Abstract type for closed convex polyhedral sets.

Notes

See HPolyhedron for a standard implementation of this interface.

Every concrete AbstractPolyhedron must define the following functions:

• constraints_list(::AbstractPolyhedron) – 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

constraints_list(X::LazySet)

Compute a list of linear constraints of a polyhedral set.

Input

• X – polyhedral set

Output

A list of the linear constraints of X.

an_element(P::AbstractPolyhedron; [solver]=default_lp_solver(eltype(P)))

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.

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.

isbounded(P::AbstractPolyhedron; [solver]=default_lp_solver(eltype(P)))

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.

isuniversal(X::LazySet, witness::Bool=false)

Check whether a set is universal.

Input

• X – set
• witness – (optional, default: false) compute a witness if activated

Output

• If the witness option is deactivated: true iff $X = ℝ^n$
• If the witness option is activated:
  • (true, []) iff $X = ℝ^n$
  • (false, v) iff $X ≠ ℝ^n$ for some $v ∉ X$

Extended help

isuniversal(P::AbstractPolyhedron, [witness]::Bool=false)

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.

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

∈(x::AbstractVector, X::LazySet)

Check whether a point lies in a set.

Input

• x – point/vector
• X – set

Output

true iff $x ∈ X$.

Extended help

∈(x::AbstractVector, P::AbstractPolyhedron)

Algorithm

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

project(P::AbstractPolyhedron, block::AbstractVector{Int}; [kwargs...])

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.
• 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)

intersection(P1::AbstractPolyhedron{N}, P2::AbstractPolyhedron{N};
             [backend]=default_lp_solver(N), [prune]::Bool=true) where {N}

Compute the intersection of two polyhedra.

Input

• P1 – polyhedron
• P2 – polyhedron
• backend – (optional, default: default_lp_solver(N)) the LP solver used for the removal of redundant constraints; see the Notes section below for details
• prune – (optional, default: true) flag for removing redundant constraints

Output

An HPolyhedron resulting from the intersection of P1 and P2, with the redundant constraints removed, or an empty set if the intersection is empty. If one of the arguments is a polytope, the result is an HPolytope instead.

Notes

The default value of the solver backend is default_lp_solver(N) and it is used to run a feasiblity LP to remove the redundant constraints of the intersection.

If you want to use the Polyhedra library, pass an appropriate backend. For example, use default_polyhedra_backend(P) for the default Polyhedra library, or use CDDLib.Library() for the CDD library.

There are some shortcomings of the removal of constraints using the default Polyhedra library; see e.g. #1038 and Polyhedra#146. It is safer to check for emptiness of intersection before calling this function in those cases.

Algorithm

This implementation unifies the constraints of the two sets obtained from the constraints_list method.

minkowski_sum(P::AbstractPolyhedron, Q::AbstractPolyhedron;
              [backend]=nothing, [algorithm]=nothing, [prune]=true)

Compute the Minkowski sum of two polyhedra in constraint representation.

Input

• P – polyhedron in constraint representation
• Q – polyhedron in constraint representation
• backend – (optional, default: nothing) polyhedral computations backend
• algorithm – (optional, default: nothing) algorithm to eliminate variables; available options are Polyhedra.FourierMotzkin, Polyhedra.BlockElimination, and Polyhedra.ProjectGenerators
• prune – (optional, default: true) if true, apply a post-processing to remove redundant constraints

Output

A polyhedron in H-representation that corresponds to the Minkowski sum of P and Q.

Algorithm

This function implements the concrete Minkowski sum by projection and variable elimination as detailed in Kvasnica [Kva05]. The idea is that if we write $P$ and $Q$ in simple H-representation, that is, $P = \{x ∈ ℝ^n : Ax ≤ b \}$ and $Q = \{x ∈ ℝ^n : Cx ≤ d \}$, then their Minkowski sum can be seen as the projection onto the first $n$-dimensional coordinates of the polyhedron:

\[ \begin{pmatrix} 0 & A \ C & -C \end{pmatrix} \binom{x}{y} ≤ \binom{b}{d}\]

This is seen by noting that $P ⊕ Q$ corresponds to the set of points $x ∈ ℝ^n$ such that $x = y + z$ with $Ay ≤ b$ and $Cz ≤ d$; hence it follows that $Ay ≤ b$ and $C(x-y) ≤ d$, and the inequality above follows by considering the $2n$-dimensional space $\binom{x}{y}$. The reduction from $2n$ to $n$ variables is performed using an elimination algorithm as described next.

The elimination of variables depends on the polyhedra library Polyhedra, which itself uses CDDLib for variable elimination. The available algorithms are:

• Polyhedra.FourierMotzkin – projection by computing the H-representation and applying the Fourier-Motzkin elimination algorithm to it

• Polyhedra.BlockElimination – projection by computing the H-representation and applying the block elimination algorithm to it

• Polyhedra.ProjectGenerators – projection by computing the V-representation

_isbounded_stiemke(constraints::AbstractVector{<:HalfSpace{N}};
                   solver=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., Mangasarian [Man94].

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

_linear_map_polyhedron(M::AbstractMatrix,
                       P::LazySet;
                       [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)

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

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.

addconstraint!(P::AbstractPolyhedron, constraint::HalfSpace)

Add a linear constraint to a set in constraint representation in-place.

Input

• P – set in constraint representation
• constraint – linear constraint to add

Notes

It is left to the user to guarantee that the dimension of all linear constraints is the same.

ishyperplanar(P::AbstractPolyhedron)

Determine whether a polyhedron is equivalent to a hyperplane.

Input

• P – polyhedron

Output

true iff P is hyperplanar, i.e., consists of two linear constraints $a·x ≤ b$ and $-a·x ≤ -b$.

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.