Optimizer returns status infeasible for a feasible problem

I’m running into the following issue while solving a set of LP-problems using JuMP:

I want to create a model, optimize it, and then in turn minimize and maximize its variables while requiring the original optimum to be achieved. For those in the know, yes, it’s FVA. The code would look something like this

using JuMP, CPLEX
model = JuMP.Model(CPLEX.Optimizer)
S = rand(5,5)
lower = zeros(5)
upper = ones(5)
c = [1; ones(4)]'
b = zeros(5)
@variable(model, x[i=1:5], lower_bound=lower[i], upper_bound=upper[i])
@objective(model, Max, c * x)
@constraint(model, S * x .== b)
optimum = JuMP.objective_value(model)
@constraint(model, c * x ≥ optimum)

for i in 1:5
      @objective(model, MOI.MIN_SENSE, x[i])
      min = JuMP.objective_value(model)

      JuMP.set_objective_sense(model, MOI.MAX_SENSE)
      max = JuMP.objective_value(model)

Now by definition, if the original problem is feasible, all of the subproblems within the for-loop should be as well. However, it appears that sometimes the optimizer enters into a “degenerate” state, and return infeasible for one of the subproblems. If a new optimizer instance is given, it returns feasible as should. In other words

# termination_status(model) will return
# INFEASIBLE::TerminationStatusCode = 2

set_optimizer(model, CPLEX.Optimizer);
# termination_status(model) will return
# OPTIMAL::TerminationStatusCode = 1

I’m at a loss as to what is going on, I’m not even sure at what level the problem is. I originally thought the problem was related to CPLEX but it appears it happens with Gurobi as well (but not necessarily with the exact same model).

I’m on
Julia Version 1.3.0
Commit 46ce4d7933 (2019-11-26 06:09 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)

JuMP v0.21.3
CPLEX v0.6.6

Is that a typo in termination_status(optimization_model)? It should say model, right? I just want to make sure.

1 Like

Please read the first post of Please read: make it easier to help you. You should format your code, and provide a minimal working example. I can’t run your code because I don’t know what lower, upper, c, b, and S are.

I’ve seen infeasibilities like this occur because the presolve tolerance of the solver is different to the final tolerance. (So using an optimal solution in the next solve may be classified as infeasible.)

I suggest you try @constraint(model, c * x ≥ optimum - 1e-4) (or some other similar tolerance).