Concrete Polyhedra

The focus of LazySets.jl is to wrap set representations and operations into specialized types, delaying the evaluation of the result of an expression until it is necessary. However, sometimes it is desirable to perform an explicit computation. For concrete operations with polyhedra we rely on the polyhedra manipulation library Polyhedra.jl.

Polyhedra.jl provides a unified interface to well-known implementations of polyhedral computations (which we also call backends) such as CDD or LRS. See the complete list in the documentation of Polyhedra.jl. This is a great advantage because we can easily use a library that supports floating point arithmetic, rational arithmetic, multiple precision, etc. The libraries also include projection and elimination of variables through Fourier-Motzkin. If you are interested in specific numeric types different from the default Float64, such as Float32, these may not be supported by the backend, in which case Julia may automatically promote to, e.g., Float64. As an example, CDD (which is used via the wrapper package CDDLib.jl) can only be used with numeric type Float64 for floating-point arithmetic and with numeric type Rational{BigInt} for exact arithmetic.

Below we give examples of operations that are performed using Polyhedra.jl.

• Concrete Polyhedra
• • Creating polyhedra
  • Methods
  • Projections

Creating polyhedra

To use the Polyhedra.jl interface, you need to load the package with using Polyhedra. Let's create an H-representation object:

using Plots, LazySets, Polyhedra, LinearAlgebra

A = [1. 1;1 -1;-1 0]
b = [1.,0,0]
H = Polyhedra.hrep(A, b)

H-representation Polyhedra.MixedMatHRep{Float64, Matrix{Float64}}:
3-element iterator of Polyhedra.HalfSpace{Float64, Vector{Float64}}:
 HalfSpace([1.0, 1.0], 1.0)
 HalfSpace([1.0, -1.0], 0.0)
 HalfSpace([-1.0, 0.0], 0.0)

It is used to instantiate a new polyhedron:

p = polyhedron(H)

Polyhedron Polyhedra.DefaultPolyhedron{Float64, Polyhedra.MixedMatHRep{Float64, Matrix{Float64}}, Polyhedra.MixedMatVRep{Float64, Matrix{Float64}}}:
3-element iterator of Polyhedra.HalfSpace{Float64, Vector{Float64}}:
 HalfSpace([1.0, 1.0], 1.0)
 HalfSpace([1.0, -1.0], 0.0)
 HalfSpace([-1.0, 0.0], 0.0)

Now, p is of the generic type Polyhedra.SimplePolyhedron{2,Float64, ...}, where 2 states for its ambient dimension, and Float64 the numeric field. The remaining fields specify the type of representation:

typeof(p)

Polyhedra.DefaultPolyhedron{Float64, Polyhedra.MixedMatHRep{Float64, Matrix{Float64}}, Polyhedra.MixedMatVRep{Float64, Matrix{Float64}}}

Observe that we can use a particular backend, such as the CDD library:

using CDDLib

p = polyhedron(H, CDDLib.Library())

Polyhedron CDDLib.Polyhedron{Float64}:
3-element iterator of Polyhedra.HalfSpace{Float64, Vector{Float64}}:
 HalfSpace([1.0, 1.0], 1.0)
 HalfSpace([1.0, -1.0], 0.0)
 HalfSpace([-1.0, 0.0], 0.0)

On the other hand, a LazySets.HPolytope object can be constructed from p:

x = HPolytope(p)
x.constraints

3-element Vector{HalfSpace{Float64, Vector{Float64}}}:
 HalfSpace{Float64, Vector{Float64}}([1.0, 1.0], 1.0)
 HalfSpace{Float64, Vector{Float64}}([1.0, -1.0], 0.0)
 HalfSpace{Float64, Vector{Float64}}([-1.0, 0.0], 0.0)

Conversely, from a HPolytope we can build a polyhedron:

y = polyhedron(x)
typeof(y)

Polyhedra.DefaultPolyhedron{Float64, Polyhedra.MixedMatHRep{Float64, Matrix{Float64}}, Polyhedra.MixedMatVRep{Float64, Matrix{Float64}}}

Moreover, you can specify the backend with an extra argument. For instance, we can use an exact representation through the Library(:exact):

A, b = Rational{Int}[1 1;1 -1;-1 0], Rational{Int}[1,0,0]
p = HPolytope(A, b)

polyhedron(p; backend=CDDLib.Library(:exact))

Polyhedron CDDLib.Polyhedron{Rational{BigInt}}:
3-element iterator of Polyhedra.HalfSpace{Rational{BigInt}, Vector{Rational{BigInt}}}:
 HalfSpace(Rational{BigInt}[1, 1], 1//1)
 HalfSpace(Rational{BigInt}[1, -1], 0//1)
 HalfSpace(Rational{BigInt}[-1, 0], 0//1)

Methods

The utility methods available are convex hull, intersection and cartesian product. The dual representation as a list of vertices can be obtained with the vertices_list function.

p = HPolytope([LinearConstraint([1.0, 0.0], 1.0),
               LinearConstraint([0.0, 1.0], 1.0),
               LinearConstraint([-1.0, 0.0], 1.0),
               LinearConstraint([0.0, -1.0], 1.0)])

constraints_list(p)

4-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)
 HalfSpace{Float64, Vector{Float64}}([-1.0, 0.0], 1.0)
 HalfSpace{Float64, Vector{Float64}}([0.0, -1.0], 1.0)

vertices_list(p)

4-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [-1.0, 1.0]
 [-1.0, -1.0]
 [1.0, -1.0]

For example, the concrete intersection of two polytopes is performed with the intersection method.

E = Ellipsoid(ones(2), Diagonal([2.0, 0.5]))
B = Ball1([2.5, 1.5], .8)

import LazySets.Approximations.overapproximate
polyoverapprox(x) = HPolytope(overapproximate(x, 1e-3).constraints)

Epoly = polyoverapprox(E)
Bpoly = polyoverapprox(B)
X = intersection(Epoly, Bpoly)

plot(E, 1e-3, aspectratio=1, alpha=0.4)
plot!(B, 1e-3, alpha=0.4)
plot!(X, 1e-3, alpha=0.4, color="black")

Projections

Projection of high-dimensional polyhedra and elimination of variables can be performed with the eliminate function, which supports three types of methods: :FourierMotzkin, :BlockElimination and :ProjectGenerators.

For further details, see the documentation of Polyhedra.jl.