Objective. This notebook contains illustrating examples for the utilities in the polyhedron_tools module.
In [1]:
%display typeset
In [2]:
from polyhedron_tools.misc import polyhedron_from_Hrep, polyhedron_to_Hrep
A = matrix([[-1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[ 0.0, 1.0, 0.0, 0.0, 0.0, 0.0],
[ 0.0, -1.0, 0.0, 0.0, 0.0, 0.0],
[ 0.0, 0.0, -1.0, 0.0, 0.0, 0.0],
[ 0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
[ 0.0, 0.0, 0.0, -1.0, 0.0, 0.0],
[ 0.0, 0.0, 0.0, 1.0, 0.0, 0.0],
[ 0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
[ 0.0, 0.0, 0.0, 0.0, -1.0, 0.0],
[ 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
[ 0.0, 0.0, 0.0, 0.0, 0.0, -1.0]])
b = vector([0.0, 10.0, 0.0, 0.0, 0.2, 0.2, 0.1, 0.1, 0.0, 0.0, 0.0, 0.0])
P = polyhedron_from_Hrep(A, b);
P
Out[2]:
In [3]:
P.inequalities()
Out[3]:
In [4]:
P.equations()
Out[4]:
It is possible to obtain the matrices that represent the inequality and the equality constraints separately, using the keyword argument separate_equality_constraints. This type of information is somtimes useful for optimization solvers.
In [5]:
[A, b, Aeq, beq] = polyhedron_to_Hrep(P, separate_equality_constraints = True)
In [6]:
A, b
Out[6]:
In [7]:
Aeq, beq
Out[7]:
Let's construct a ball in the infinity norm, specifying the center and radius.
We remark that the case of a hyperbox can be done in Sage's library polytopes.hypercube(n) where n is the dimension. However, as of Sage v7.6. there is no such hyperrectangle function (or n-orthotope, see the wikipedia page), so we use BoxInfty.
In [8]:
from polyhedron_tools.misc import BoxInfty
In [9]:
P = BoxInfty(center=[1,2,3], radius=0.1); P.plot(aspect_ratio=1)
Out[9]:
As a side note, the function also works when the arguments are not named, as in
In [10]:
P = BoxInfty([1,2,3], 0.1); P
Out[10]:
Another use of BoxInfty is to specify the lengths of the sides. For example:
In [12]:
P = BoxInfty([[0,1], [2,3]]); P
Out[12]:
The random_polygon_2d receives the number of arguments as input and produces a polygon whose vertices are randomly sampled from the unit circle. See the docstring for more another options.
In [13]:
from polyhedron_tools.polygons import random_polygon_2d
random_polygon_2d(5)
Out[13]:
In [14]:
from polyhedron_tools.misc import BoxInfty, opposite_polyhedron
P = BoxInfty([1,1], 0.5);
mp = opposite_polyhedron(P);
In [15]:
P.plot(aspect_ratio=1) + mp.plot(color='red')
Out[15]:
In [16]:
from polyhedron_tools.misc import support_function
P = BoxInfty([1,2,3,4,5], 1); P
support_function(P, [1,-1,1,-1,1], verbose=1)
Out[16]:
It is also possible to input the polyhedron in matrix form, $[A, b]$. If this is possible, it is preferable, since it is often faster. Below is an example with $12$ variables. We get beteen 3x and 4x improvement in the second case.
In [17]:
reset('P, A, b')
P = BoxInfty([1,2,3,4,5,6,7,8,9,10,11,12], 1); P
[A, b] = polyhedron_to_Hrep(P)
In [18]:
timeit('support_function(P, [1,-1,1,-1,1,-1,1,-1,1,-1,1,-1])')
In [19]:
support_function([A, b], [1,-1,1,-1,1,-1,1,-1,1,-1,1,-1])
Out[19]:
In [20]:
timeit('support_function([A, b], [1,-1,1,-1,1,-1,1,-1,1,-1,1,-1])')
In [21]:
from polyhedron_tools.misc import support_function, support_function_ellipsoid
import random
# Generate a random ellipsoid and check support function outer approximation.
# Define an ellipse as: x^T*Q*x <= 1
M = random_matrix(RR, 2, distribution="uniform")
Q = M.T*M
f = lambda x, y : Q[0,0]*x^2 + Q[1,1]*y^2 + (Q[0,1]+Q[1,0])*x*y-1
E = implicit_plot(f,(-5,5),(-3,3),fill=True,alpha=0.5,plot_points=600)
# generate at random k directions, and compute the overapproximation of E using support functions
# It works 'in average': we might get unbounded domains (random choice did not enclose the ellipsoid).
# It is recommended to use QQ as base_ring to avoid 'frozen set' issues.
k=15
A = matrix(RR,k,2); b = vector(RR,k)
for i in range(k):
theta = random.uniform(0, 2*pi.n(digits=5))
d = vector(RR,[cos(theta), sin(theta)])
s_fun = support_function_ellipsoid(Q, d)
A.set_row(i,d); b[i] = s_fun
OmegaApprox = polyhedron_from_Hrep(A, b, base_ring = QQ)
E + OmegaApprox.plot(fill=False, color='red')
Out[21]:
In [22]:
from polyhedron_tools.misc import BoxInfty, radius
In [23]:
P = BoxInfty([-13,24,-51,18.54,309],27.04);
radius(P)
In [24]:
got_lengths, got_center_and_radius = False, False
In [25]:
got_lengths is not False
Out[25]:
In [26]:
radius(polyhedron_to_Hrep(P))
Out[26]:
In [28]:
8401/25.n()
Out[28]:
Consider a higher-dimensional system. We obtain almost a 200x improvement for a 15-dimensional set. This is because in the case we call poly_sup_norm with a polytope, the bounding_box() function consumes time.
In [30]:
%%time
P = BoxInfty([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1], 14.28);
radius(P)
In [31]:
%%time
[A, b] = BoxInfty([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1], 14.28, base_ring = RDF, return_HSpaceRep = True)
radius([A, b])
The improvement in speed is quite interesting! Moreover, for a 20-dimensional set, the polyhedron construct does not finish.
The operation matrix x polyhedron can be performed in Sage with the usual command, *.
In [33]:
U = BoxInfty([[0,1],[0,1]])
B = matrix([[1,1],[1,-1]])
P = B * U
P.plot(color='red', alpha=0.5) + U.plot()
Out[33]:
In [36]:
from polyhedron_tools.misc import chebyshev_center, BoxInfty
from polyhedron_tools.polygons import random_polygon_2d
P = random_polygon_2d(10, base_ring = QQ)
c = chebyshev_center(P);
B = BoxInfty([[1,2],[0,1]])
b = chebyshev_center(B)
fig = point(c, color='blue') + P.plot(color='blue', alpha=0.2)
fig += point(b, color='red') + B.plot(color='red', alpha=0.2)
fig += point(P.center().n(), color='green',marker='x')
fig += point(B.center().n(), color='green',marker='x')
fig
Out[36]:
The method center() existent in the Polyhedra class, computes the average of the vertices. In contrast, the Chebyshev center is the center of the largest box enclosed by the polytope.
In [37]:
B.bounding_box()
Out[37]:
In [38]:
P.bounding_box()
Out[38]:
In [39]:
e = [ (P.bounding_box()[0][0] + P.bounding_box()[1][0])/2, (P.bounding_box()[0][1] + P.bounding_box()[1][1])/2]
l = [[P.bounding_box()[0][0], P.bounding_box()[1][0]], [P.bounding_box()[0][1], P.bounding_box()[1][1]] ]
fig += point(e,color='black') + BoxInfty(lengths=l).plot(alpha=0.1,color='grey')
fig
Out[39]:
Here we have added in grey the bounding box that is obtained from the method bounding_box(). To make the picture complete, we should also add the box of center Cheby center, and of maximal radius which is included in it.
In [41]:
from polyhedron_tools.projections import lotov_algo
from polyhedron_tools.misc import polyhedron_to_Hrep, polyhedron_from_Hrep
A = matrix([[-1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[ 0.0, 1.0, 0.0, 0.0, 0.0, 0.0],
[ 0.0, -1.0, 0.0, 0.0, 0.0, 0.0],
[ 0.0, 0.0, -1.0, 0.0, 0.0, 0.0],
[ 0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
[ 0.0, 0.0, 0.0, -1.0, 0.0, 0.0],
[ 0.0, 0.0, 0.0, 1.0, 0.0, 0.0],
[ 0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
[ 0.0, 0.0, 0.0, 0.0, -1.0, 0.0],
[ 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
[ 0.0, 0.0, 0.0, 0.0, 0.0, -1.0]])
b = vector([0.0, 10.0, 0.0, 0.0, 0.2, 0.2, 0.1, 0.1, 0.0, 0.0, 0.0, 0.0])
P = polyhedron_from_Hrep(A, b); P
Out[41]:
In [42]:
lotov_algo(A, b, [1,0,0], [0,1,0], 0.5)
Out[42]: