Best way of checking whether a polyhedron contains an integral point

Specific Domains Optimization (Mathematical)
,
1

Hey!

I’m quite new to Julia and try to solve the following problem efficiently:

For a given bounded polyhedron (by the linear equations Ax = 0 and bounds on x), I have to check whether it contains a non-negative integral point except the zero vector (which is always feasible). So far, I’m using JuMP to maximize something like sum(x) and check if the optimal value is 0 or greater than 0.

Now I’ve found the Polyhedra package, but no function that does what I need… Has anyone an idea how this can be implemented more efficiently than I’ve described above?

Thank you very much for any hint!

2

Your approach with JuMP sounds OK to me.

You can also omit the objective function (just solve find any feasible point) and avoid the origin with a constraint of the form sum(x) >= 1, right?

3

Thank you for your answer.

Yeah right, I can omit the objective function and add a constraint to avoid the origin, thanks!

Adittional question: For some instances I obtain the termination_status DUAL_INFEASIBLE. Is there any chance to decided feasibility of the primal anyway?

4

What solver are you using? DUAL_INFEASIBLE suggest that the problem is unbounded, but that shouldn’t matter for a feasibility problem? Maybe you can provide a bounding box as well?

5

I’m using CPLEX and Gurobi with the same result (DUAL_INFEASIBLE). That was also my first thought but the problem is bounded… I can offer the example below (sorry, I haven’t found a smaller matrix where this phenomenon occur).

If I switch to GLPK, I obtain warnings like Warning: numerical instability (primal simplex, phase I) and Warning: numerical instability (primal simplex, phase II), but the termination status is now OPTIMAL in this specific example… I don’t really get it.

using CPLEX
using JuMP

A = [1 1 1 1 -1 0 0 0 -1 0 0 0 0 0 -1 0 0 0 0 0 -1 0 0 0 0 0 0 0;
     -1 0 0 0 1 1 1 1 0 0 0 -1 0 0 0 0 0 -1 0 0 0 0 0 0 -1 0 0 0;
     0 -1 0 0 0 0 0 0 1 1 1 0 0 0 0 -1 0 0 0 0 0 -1 0 0 0 0 0 0;
     0 0 0 0 0 -1 0 0 0 0 0 1 1 1 0 0 0 0 -1 0 0 0 0 0 0 -1 0 0;
     0 0 -1 0 0 0 0 0 0 -1 0 0 0 0 1 1 1 0 0 0 0 0 -1 0 0 0 0 0;
     0 0 0 0 0 0 -1 0 0 0 0 0 -1 0 0 0 0 1 1 1 0 0 0 0 0 0 -1 0;
     0 0 0 -1 0 0 0 0 0 0 -1 0 0 0 0 0 -1 0 0 0 1 1 1 1 0 0 0 -1;
     0 0 0 0 0 0 0 -1 0 0 0 0 0 -1 0 0 0 0 0 -1 0 0 0 -1 1 1 1 1;
     1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1;
     0 -1 0 0 1 1 1 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
     0 0 -1 0 0 -1 0 0 1 1 1 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0;
     0 0 0 -1 0 0 -1 0 0 -1 0 1 1 1 0 -1 0 -1 0 0 -1 0 0 0 0 0 0 0;
     0 0 0 0 0 0 0 -1 0 0 -1 0 -1 0 1 1 1 0 -1 0 0 -1 0 0 -1 0 0 0;
     0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 1 1 1 0 0 -1 0 0 -1 0 0;
     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 1 1 1 0 0 -1 0;
     -1 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1]

nrows, ncols = size(A,1), size(A,2)

M = Model(with_optimizer(CPLEX.Optimizer, CPX_PARAM_SCRIND=false))

@variable(M, x[1:ncols], lower_bound=0, upper_bound=(ncols + 1)*ncols^(ncols/2), Int)
@objective(M, Max, sum(x))
@constraint(M, con, A*x .== zeros(nrows))

optimize!(M)
println(termination_status(M))
6

That bound is rather wide. It evaluates to > 5e21 for me. I’m guessing that most solvers consider this value “infinite”, hence the variables are actually unbounded.

Also, your example still has the @objective set.

When reading the first post, I was thinking of binary variables, for some reason.

7

Ah, that makes sense. I changed the example as you suggested, and it works now! (Probably because the problem cannot be unbounded anymore?)

Thank you very much for your help!!
(I’m not sure which answer I should mark as solution to my question…)

8

When calling isempty on a polyhedron with Polyhedra, it solves the same linear problem with JuMP than the one you are formulating in JuMP so you won’t gain anything by using Polyhedra, it’s the same for the solvers

2 Likes
9

Thanks for your answer! That is exactly the answer I was looking for, i.e., if Polyhedra can do it better than me…