# Best way of checking whether a polyhedron contains an integral point

**URL:** <https://discourse.julialang.org/t/best-way-of-checking-whether-a-polyhedron-contains-an-integral-point/36766>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump\
**Created:** [March 30, 2020, 10:44pm UTC](https://discourse.julialang.org/t/best-way-of-checking-whether-a-polyhedron-contains-an-integral-point/36766 "2020-03-30T22:44:14Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![python15](https://avatars.discourse-cdn.com/v4/letter/p/ad7895/32.png) [@python15](https://discourse.julialang.org/u/python15)\
**Post date:** [March 30, 2020, 10:44pm UTC](https://discourse.julialang.org/t/best-way-of-checking-whether-a-polyhedron-contains-an-integral-point/36766/1 "2020-03-30T22:44:14Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![leethargo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leethargo/32/6004_2.png) [@leethargo](https://discourse.julialang.org/u/leethargo)\
**Post date:** [March 31, 2020, 5:44am UTC](https://discourse.julialang.org/t/best-way-of-checking-whether-a-polyhedron-contains-an-integral-point/36766/2 "2020-03-31T05:44:33Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![python15](https://avatars.discourse-cdn.com/v4/letter/p/ad7895/32.png) [@python15](https://discourse.julialang.org/u/python15)\
**Post date:** [March 31, 2020, 9:09am UTC](https://discourse.julialang.org/t/best-way-of-checking-whether-a-polyhedron-contains-an-integral-point/36766/3 "2020-03-31T09:09:13Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![leethargo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leethargo/32/6004_2.png) [@leethargo](https://discourse.julialang.org/u/leethargo)\
**Post date:** [March 31, 2020, 10:13am UTC](https://discourse.julialang.org/t/best-way-of-checking-whether-a-polyhedron-contains-an-integral-point/36766/4 "2020-03-31T10:13:34Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![python15](https://avatars.discourse-cdn.com/v4/letter/p/ad7895/32.png) [@python15](https://discourse.julialang.org/u/python15)\
**Post date:** [March 31, 2020, 10:51am UTC](https://discourse.julialang.org/t/best-way-of-checking-whether-a-polyhedron-contains-an-integral-point/36766/5 "2020-03-31T10:51:16Z")

</div>

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.

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

```

---

<div class="post-metadata">

**Author:** ![leethargo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leethargo/32/6004_2.png) [@leethargo](https://discourse.julialang.org/u/leethargo)\
**Post date:** [March 31, 2020, 11:06am UTC](https://discourse.julialang.org/t/best-way-of-checking-whether-a-polyhedron-contains-an-integral-point/36766/6 "2020-03-31T11:06:33Z")

</div>

> [@python15](#):
>
> `upper_bound=(ncols + 1)*ncols^(ncols/2)`

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.

---

<div class="post-metadata">

**Author:** ![python15](https://avatars.discourse-cdn.com/v4/letter/p/ad7895/32.png) [@python15](https://discourse.julialang.org/u/python15)\
**Post date:** [March 31, 2020, 11:27am UTC](https://discourse.julialang.org/t/best-way-of-checking-whether-a-polyhedron-contains-an-integral-point/36766/7 "2020-03-31T11:27:22Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![blegat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blegat/32/217090_2.png) [@blegat](https://discourse.julialang.org/u/blegat)\
**Post date:** [April 2, 2020, 6:04pm UTC](https://discourse.julialang.org/t/best-way-of-checking-whether-a-polyhedron-contains-an-integral-point/36766/8 "2020-04-02T18:04:50Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![python15](https://avatars.discourse-cdn.com/v4/letter/p/ad7895/32.png) [@python15](https://discourse.julialang.org/u/python15)\
**Post date:** [April 2, 2020, 6:57pm UTC](https://discourse.julialang.org/t/best-way-of-checking-whether-a-polyhedron-contains-an-integral-point/36766/9 "2020-04-02T18:57:41Z")

</div>

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