Gurobi's optimal solution is actually infeasible

[XG Yang]

I am trying to solve a simple MILP as follows:

SETS:
i
j

PARAMETERS
A[i,j]

REAL VARIABLES:
x[j]

BINARY VARIABLES
y[j]

OBJECTIVE FUNCTION:
z = minimize(sum([y[j] for j in J]))

CONSTRAINTS
Constraint1: sum([A[i,j]*x[j] for j in J]) == 0,   for i in I
Constraint2: x[j] >= -10*y[j],       for j in J
Constraint3: x['specific j'] >= 10,    'specific j' in J

Now, when I solve this with gurobi it returns an optimal solution with an objective value of zero (or close to zero like 1e-10), i.e., at optimal solution all y[j]'s are zero. However, when i fix all y[j]'s to zero and resolve the same problem it becomes infeasible. The constraint that is violated is Constraint 3. By manual inspection of the problem it turns out that it indeed is infeasible if all y[j]'s are zero. So, I was wondering why gurobi gives this as an optimal solution and if there is any way to fix it. I set the solver parameters FeasibilityTol, OptimalityTol, MIPGap, MIPGapAbs and IntFeasTol to 1e-9. The same problem is solved perfectly fine when using different 'specific j' in Constraint 3 (i.e., this happens only for one 'specific j'

[Tobias Achterberg]

Can you send the Gurobi log output? I suspect that you have weird numbers in
your A matrix that make the model numerically instable.

Tobias

[XG Yang]

Sure. Would you let me know where the gurobi.log file is saved and how I should tell gurobi to save the model and other info in the log file? I have a gurobi.log file in the directory where my code is located and all what I see in that file is:

Gurobi 6.0.0 (linux64) logging started Tue Oct  6 17:32:33 2015

I am accessing gurobi through pyomo.

[Sonja Mars]

Hi,

How are you exactly calling Pyomo? Do you use a Python script or do you use pyomo solve?

For example, if you use a Python script you will have to add "tee=True" to the solve command. So the solve command would look for example like this:

results = opt.solve(m, tee=True)

This will produce Gurobi output on the terminal. Please also take a look here: https://software.sandia.gov/downloads/pub/pyomo/PyomoOnlineDocs.html#_display_of_solver_output

Thanks and best regards,
Sonja

-----------------------------------------------------------------
Dr. Sonja Mars
Gurobi Optimization

> --
>
> ---
> You received this message because you are subscribed to the Google Groups "Gurobi Optimization" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to gurobi+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[XG Yang]

Hello - I am using "results = opt.solve(m, tee=True)" and I know that tee=True will show the solver output but I thought Tobias is asking for a logfile wherein details of the optimization models can be found (to check for weird numbers in matrix A). The solver output for my problem is as follows:

Changed value of parameter FeasibilityTol to 1e-09
   Prev: 1e-06   Min: 1e-09   Max: 0.01   Default: 1e-06
Parameter MIPGapAbs unchanged
   Value: 1e-10   Min: 0.0   Max: 1e+100   Default: 1e-10
Changed value of parameter Threads to 4
   Prev: 0   Min: 0   Max: 1024   Default: 0
Changed value of parameter MIPGap to 1e-09
   Prev: 0.0001   Min: 0.0   Max: 1e+100   Default: 0.0001
Changed value of parameter VarBranch to 3
   Prev: -1   Min: -1   Max: 3   Default: -1
Changed value of parameter NodefileStart to 1.0
   Prev: 1e+100   Min: 0.0   Max: 1e+100   Default: 1e+100
Changed value of parameter OptimalityTol to 1e-09
   Prev: 1e-06   Min: 1e-09   Max: 0.01   Default: 1e-06
Changed value of parameter IntFeasTol to 1e-09
   Prev: 1e-05   Min: 1e-09   Max: 0.1   Default: 1e-05
Optimize a model with 2106 rows, 2883 columns and 10783 nonzeros
Coefficient statistics:
  Matrix range    [2e-06, 1e+03]
  Objective range [1e+00, 1e+00]
  Bounds range    [1e+00, 1e+03]
  RHS range       [2e-03, 1e+00]
Presolve removed 1426 rows and 1542 columns
Presolve time: 0.04s
Presolved: 680 rows, 1341 columns, 5688 nonzeros
Variable types: 1123 continuous, 218 integer (218 binary)

Root relaxation: objective 0.000000e+00, 957 iterations, 0.03 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0       0.0000000    0.00000  0.00%     -    0s

Explored 0 nodes (961 simplex iterations) in 0.09 seconds
Thread count was 4 (of 12 available processors)

Optimal solution found (tolerance 1.00e-09)
Best objective 5.001140608609e-10, best bound 5.001140608609e-10, gap 0.0%
solution =  []

Objective value = 5.00114060861e-10, Optimality status = optimal, Solution status = optimal, Solver run status = normal

I don't see anything unusual in this output.
An an update, I have verified (by using a loop) that the problem has two alterantive solutions with an objective value of one (i.e., when only one binary variables is one). The solver is able to find one of these solutions if I fix a large number of binary variables to zero.

One thing that is unclear to me is that if the problem is a due to numeric issues, then shouldn't it always behave like this? Like I mentioned in my first post, I have this issues for only one "special j" and not for others.

[Sonja Mars]

Hi,

You are right, it is not an unusual output. However, your model has numerical issues. If you take a look at the coefficients:

> Coefficient statistics:
> Matrix range [2e-06, 1e+03]
> Objective range [1e+00, 1e+00]
> Bounds range [1e+00, 1e+03]
> RHS range [2e-03, 1e+00]

You can see a large matrix range and this can cause numerical issues, and it can cause the behavior you are seeing.

Would it be possible for you to send the model file (we would prefer an lp or mps file) so that we can investigate? If you don't want to post it here, please send it to support[at]gurobi.com

[XG Yang]

Sure. I'll be happy to send you the model. I just don't know how to tell gurobi to create lp or mps files.

[Sonja Mars]

Hi,

You can simply use the write method of the model object included in pyomo. So for example, this should work if your model object is named "model".

model.write("myModel.mps")

Best regards,