find one interior point of a polyhedron Hi,

I have a bunch of inequations and I would like to know if there is a solution. What is the simplest way to achieve this in Sage ? I tried using MILP with success... but the point returned is not very fancy (often on the boundary). In an ideal world, I would like to optimize some quadratic function in order for the solution to be the nicest possible.

Here is a sample and simple example with only one inequality (other constraints are equalities) where I reproduced the output of MILP:

Constraints:
2.0 <= x_0 + x_7 + x_8 <= 2.0
2.0 <= x_1 + x_2 <= 2.0
2.0 <= x_3 + x_5 + x_9 <= 2.0
2.0 <= x_4 + x_6 <= 2.0
1.0 <= x_0 + x_2 + x_9 <= 1.0
1.0 <= x_5 + x_6 + x_8 <= 1.0
0.0 <= x_2 + x_6 <= 1.0
Variables:
x_0 is a continuous variable (min=0.0, max=2.0)
x_1 is a continuous variable (min=0.0, max=2.0)
x_2 is a continuous variable (min=0.0, max=2.0)
x_3 is a continuous variable (min=0.0, max=2.0)
x_4 is a continuous variable (min=0.0, max=2.0)
x_5 is a continuous variable (min=0.0, max=2.0)
x_6 is a continuous variable (min=0.0, max=2.0)
x_7 is a continuous variable (min=0.0, max=2.0)
x_8 is a continuous variable (min=0.0, max=2.0)
x_9 is a continuous variable (min=0.0, max=2.0)

There are plenty of nice solutions, one is

x0 = x2 = x5 = x6 = x8 = 1/3
x1 = x4 = 5/3
x3 = x7 = 4/3

But with the solver I got:

sage: p.solve()
sage: p.get_values(p,p,p,p,p,p,p,p,p)
[1, 2, 0, 2, 2, 0, 0, 0, 1]

Thanks

edit retag close merge delete

Sort by » oldest newest most voted

The solver aims at miximizing a linear form on the polytope you defined, hence it likes returning extreme or boundary points even if the linear form is not specified. If you want to find a (relative) interior point of your polytope, it suffice to build a non-degenerate convex combination of the extreme points of the polytope.

You can find extreme points as follows:

sage: p.polyhedron().vertices()
(A vertex at (0, 0, 1, 0, 1, 2, 1, 2, 1, 0),
A vertex at (0, 0, 0, 1, 0, 2, 1, 2, 2, 0),
A vertex at (0, 0, 1, 0, 1, 2, 1, 2, 0, 1),
A vertex at (0, 0, 0, 1, 0, 2, 1, 2, 0, 2),
A vertex at (1/2, 1/2, 0, 1/2, 0, 3/2, 3/2, 3/2, 2, 0),
A vertex at (0, 1, 0, 0, 1, 1, 2, 2, 1, 0),
A vertex at (1, 0, 1, 0, 0, 2, 1, 1, 2, 0),
A vertex at (1/2, 1/2, 0, 1/2, 0, 3/2, 3/2, 3/2, 0, 2),
A vertex at (0, 1, 0, 0, 1, 1, 2, 2, 0, 1),
A vertex at (1, 0, 1, 0, 0, 2, 1, 1, 0, 2))

For example, you can do:

sage: V = p.polyhedron().vertices()
sage: sum(v.vector() for v in V) / len(V)
(3/10, 3/10, 2/5, 3/10, 2/5, 17/10, 13/10, 17/10, 4/5, 4/5)

Equivalently:

sage: p.polyhedron().center()
(3/10, 3/10, 2/5, 3/10, 2/5, 17/10, 13/10, 17/10, 4/5, 4/5)

Or even:

sage: sum((i+1)*v.vector() for (i,v) in enumerate(V)) / (len(V)*(len(V)+1)/2)
(47/110, 43/110, 21/55, 5/22, 19/55, 177/110, 153/110, 173/110, 7/11, 56/55)
more

If all you want is some point that would be in the interior of the solution set, you can get one by doing:

sage: q = p.polyhedron()
sage: q.representative_point()

which returns the barycenter of the vertices of the polyhedron and adds the potential generating rays. So it would be in the relative interior of the polyhedron.

more