Tolerance issue for gurobi leading to incorrect/out-of-domain assignments

[hbierlee]

In the following model, we get integer variables with domains 1..3 assigned to 0 by the Gurobi interface, since Gurobi assigns them e.g. 0.999999, and we use int to convert values from Gurobi, which floors/truncates them.

import cpmpy as cp

x = cp.intvar(1, 3, shape=2, name="x")
p = cp.intvar(0, 1, shape=3, name="p")
q = cp.intvar(0, 1, shape=3, name="q")
m = cp.Model()
m += cp.sum([p[0], p[1], p[2]]) == 1
m += cp.sum([3, 3, 1, -1] * cp.cpm_array([q[0], q[1], q[2], x[1]])) == 0
m += cp.sum([q[0], q[1], q[2]]) == 1

m.solveAll(solver="gurobi", solution_limit=1000, display=x)

Outputs:

❯ python examples/tols.py 
[None 1]
[None 3]
[None 3]
[None 3]
[None 3]
[None 3]
[None 0]
[None 3]
[None 0]

Instead, we should simply use e.g. round to get more correct results.

For additional safety, we could add an integrality check (as is done for the objective, or it could be done using np.isclose as proposed in #579). Currently, we silently covert the objective to a float if it is not integer. We could:

a) do the same for integer variables, or
b) use python round in any case, or
c) round using the currently set IntFeasTol (there is no Gurobi-official way to do this, but I've seen mention of np.ceil(x - 0.5 + IntFeasTol))

Based on my (recently developed) understanding of tolerances, with e.g. the default IntFeasTol set to 1e-5, we should never see e.g. 0.98, and if we do, something is going wrong and relaxing/tightening the tolerances would not resolve it. So, in that case, for option b/c, we could additionally raise an error if the integrality check fails (although I'm not sure what advice to give to the user in this case).

PS: a different issue is that x[0] is assigned to None since it does not occur in any constraint. This can be problematic for e.g. solution counts.

[hbierlee]

Slightly updated understanding: we are guaranteed that the integer is assigned at most IntFeasTol (default 1e-5) away from an integer value. So round should always work, and a value outside the tolerance (erg. 10.1) should either raise an error (because it must be a bug), or be rounded too (risking bugs/incorrect results), or rounded but with a warning. The numpy method is for when you need to round real float decision variables.

We may still have that gurobi wrongly assigns integer variable x to 10.99999999 even if we have a constraint x<=10, as the numerical instability can shift the value of x by an arbitrary amount, tipping over into the next integer. Avoiding that requires stricter tolerances or changing constraints, or using IntegralityFocus parameter in integer problems (https://docs.gurobi.com/projects/optimizer/en/current/reference/parameters.html#integralityfocus). This can be tackled in a different issue, perhaps once we see a real case of this for our model.