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 necessary to do an explicit computation. For concrete operations with polyhedra we rely on the polyhedra manipulation library Polyhedra.jl.

Actually, Polyhedra.jl provides a unified interface to well-known implementations of polyhedral computations, 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.

Below we give examples of operations that are actually done via Polyhedra.jl.

Concrete Polyhedra

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, Compat.LinearAlgebra

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

H-representation Polyhedra.MixedMatHRep{Float64,Array{Float64,2}}:
3-element iterator of Polyhedra.HalfSpace{Float64,Array{Float64,1}}:
 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,Array{Float64,2}},Polyhedra.MixedMatVRep{Float64,Array{Float64,2}}}:
3-element iterator of Polyhedra.HalfSpace{Float64,Array{Float64,1}}:
 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,Array{Float64,2}},Polyhedra.MixedMatVRep{Float64,Array{Float64,2}}}

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,Array{Float64,1}}:
 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 Array{HalfSpace{Float64},1}:
 HalfSpace{Float64}([1.0, 1.0], 1.0)
 HalfSpace{Float64}([1.0, -1.0], 0.0)
 HalfSpace{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,Array{Float64,2}},Polyhedra.MixedMatVRep{Float64,Array{Float64,2}}}

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},Array{Rational{BigInt},1}}:
 HalfSpace(Rational{BigInt}[1//1, 1//1], 1//1)
 HalfSpace(Rational{BigInt}[1//1, -1//1], 0//1)
 HalfSpace(Rational{BigInt}[-1//1, 0//1], 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 Array{HalfSpace{Float64},1}:
 HalfSpace{Float64}([1.0, 0.0], 1.0)
 HalfSpace{Float64}([0.0, 1.0], 1.0)
 HalfSpace{Float64}([-1.0, 0.0], 1.0)
 HalfSpace{Float64}([0.0, -1.0], 1.0)

vertices_list(p)

4-element Array{Array{Float64,1},1}:
 [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.