Gurobi founds optimal solution but it is not feasible

Question

I'm currently solving a MIP model with pyomo using gurobi and I am facing strange results.

I have one constraint that looks like this:

def satisfacer_servicios_r1(modelo, s):
    return (sum(modelo.x[bv, s] for bv in modelo.BV) == 1)
modelo.satisfacer_servicios_r1 = pyo.Constraint(modelo.S, rule=satisfacer_servicios_r1)

And x is defined as:

modelo.x = pyo.Var(modelo.BV, modelo.S, domain=pyo.Binary) 

In the solution found by gurobi all x are equal to -1 nonetheless it says that constraint "satisfacer_servicios_r1" is active.

Code to check x values:

modelo.x.pprint()
x : Size=4250, Index=x_index
    Key       : Lower : Value : Upper : Fixed : Stale : Domain
       (1, 1) :     0 :  -0.0 :     1 : False : False : Binary
       (1, 5) :     0 :  -0.0 :     1 : False : False : Binary
      (1, 11) :     0 :  -0.0 :     1 : False : False : Binary
      (1, 15) :     0 :  -0.0 :     1 : False : False : Binary
      (1, 18) :     0 :  -0.0 :     1 : False : False : Binary
      (1, 31) :     0 :  -0.0 :     1 : False : False : Binary
      (1, 34) :     0 :  -0.0 :     1 : False : False : Binary

Code to check constraint status:

modelo.satisfacer_servicios_r1.pprint()
satisfacer_servicios_r1 : Size=85, Index=S, Active=True
    Key : Lower : Body                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         : Upper : Active
      1 :   1.0 :x[1,1] + x[2,1] + x[3,1] + x[4,1] + x[5,1] + x[6,1] + x[7,1] + x[8,1] + x[9,1] + x[10,1] + x[11,1] + x[12,1] + x[13,1] + x[14,1] + x[15,1] + x[16,1] + x[17,1] + x[18,1] + x[19,1] + x[20,1] + x[21,1] + x[22,1] + x[23,1] + x[24,1] + x[25,1] + x[26,1] + x[27,1] + x[28,1] + x[29,1] + x[30,1] + x[31,1] + x[32,1] + x[33,1] + x[34,1] + x[35,1] + x[36,1] + x[37,1] + x[38,1] + x[39,1] + x[40,1] + x[41,1] + x[42,1] + x[43,1] + x[44,1] + x[45,1] + x[46,1] + x[47,1] + x[48,1] + x[49,1] + x[50,1] :   1.0 :   True

Gurobi log:

x22053: 587140 rows, 21953 columns, 2317504 nonzeros
Set parameter MIPGap to value 0.15
Gurobi Optimizer version 10.0.1 build v10.0.1rc0 (win64)
CPU model: AMD Ryzen 7 5800H with Radeon Graphics, instruction set [SSE2|AVX|AVX2]
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 587140 rows, 21953 columns and 2317504 nonzeros
Model fingerprint: 0x7427a6f5
Variable types: 8653 continuous, 13300 integer (13300 binary)
Coefficient statistics:
  Matrix range     [8e-02, 5e+02]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+03]
Presolve removed 196104 rows and 8802 columns (presolve time = 5s) ...
Presolve removed 200354 rows and 8802 columns (presolve time = 10s) ...
Presolve removed 200354 rows and 8802 columns
Presolve time: 13.18s
Presolved: 386786 rows, 13151 columns, 1221701 nonzeros
Variable types: 151 continuous, 13000 integer (13000 binary)
Found heuristic solution: objective 118.1166667
Root simplex log...
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    1.0250000e+02   0.000000e+00   8.600000e+00     25s
     107    1.0250000e+02   0.000000e+00   0.000000e+00     25s
     107    1.0250000e+02   0.000000e+00   0.000000e+00     25s
     107    1.0250000e+02   0.000000e+00   0.000000e+00     25s
Use crossover to convert LP symmetric solution to basic solution...
Root crossover log...
    8514 PPushes remaining with PInf 0.0000000e+00                25s
    1154 PPushes remaining with PInf 0.0000000e+00                30s
       0 PPushes remaining with PInf 0.0000000e+00                35s
  Push phase complete: Pinf 0.0000000e+00, Dinf 2.3045463e-11     35s
Root simplex log...
Iteration    Objective       Primal Inf.    Dual Inf.      Time
    8624    1.0250000e+02   0.000000e+00   0.000000e+00     35s
Root relaxation: objective 1.025000e+02, 8624 iterations, 12.16 seconds (10.88 work units)
Total elapsed time = 60.91s
Total elapsed time = 67.15s
Total elapsed time = 73.27s
Total elapsed time = 75.11s
    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
     0     0  102.50000    0   76  118.11667  102.50000  13.2%     -   76s
Explored 1 nodes (45160 simplex iterations) in 76.75 seconds (82.07 work units)
Thread count was 16 (of 16 available processors)
Solution count 1: 118.117 
Optimal solution found (tolerance 1.50e-01)
Best objective 1.181166666666e+02, best bound 1.025000000000e+02, gap 13.2214%

Edit: I have fixed the value of one x to 1 and then check the behaviour after this. The fix works and the value is 1 on the solution, but now I have noticed that some of the x values are not 0, but -0. I do not know what this means.

Answer 1

This problem usually occurs when you have issues with numerical precision, stemming from parameters that are much larger than they should be. Things you should ask yourself:

Am I using big M? If yes, is it properly scaled? Are the parameters I'm using on very different scales? Do I have any parameters in the billions and others in other ranges? If so, it might be reasonable to consider rounding.