Best way of checking whether a polyhedron contains an integral point

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!

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?

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?

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?

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))

``````

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.

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…)

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

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