Polyhedra (AbstractPolyhedron)

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

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

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.

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.

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.

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.

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

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

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

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

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.

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.

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.

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.

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.