Conversation
From a previous solve, z_k may be a value like z+eps where z in Z and eps is the feasibility tolerance. However, if scalars[k] is integer valued, then presolve may (somewhat reasonably) deduce that this problem is infeasible. Instead of rounding z_k, changing EqualTo to LessThan seemed to work on the instances I have. I don't know why. There is also a weird situation in which Gurobi declared a problem unbounded. I can't reproduce without running the entire thing, but I think it is a mix of presolve proving that there is no finite solution (because its infeasible) and yet there being a "feasible" MIP solution in memory from the previous solve. It comes down to the weird mix of tolerances in MIP starts, presolve, simplex, and this equality constraint. This is most probably a bug in Gurobi, in that it should either report optimal or infeasible. But not unbounded. I've seem something similar in SDDP.jl.
|
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #198 +/- ##
=======================================
Coverage 99.77% 99.77%
=======================================
Files 12 12
Lines 1307 1307
=======================================
Hits 1304 1304
Misses 3 3 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
@simonbowly here's something of interest. I haven't had any luck with a minimal reproducible example, but look at the second log above. It reports Here's the model we're trying to solve: model.mps.txt Gurobi 13 says it is infeasible: Gurobi 12 says it is feasible In no case should it be unbounded. |
|
Getting somewhere: julia> using JuMP, Gurobi
julia> model = read_from_file("/tmp/model.mps")
A JuMP Model
├ solver: none
├ objective_sense: MIN_SENSE
│ └ objective_function_type: AffExpr
├ num_variables: 225
├ num_constraints: 258
│ ├ AffExpr in MOI.EqualTo{Float64}: 31
│ ├ AffExpr in MOI.LessThan{Float64}: 2
│ └ VariableRef in MOI.ZeroOne: 225
└ Names registered in the model: none
julia> set_optimizer(model, Gurobi.Optimizer)
Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID to value 722777
WLS license 722777 - registered to JuMP Development
julia> c = constraint_by_name(model, "R74")
R74 : 8 x[1,1] + 7 x[2,1] + 13 x[3,1] + 14 x[4,1] + 5 x[5,1] + 12 x[6,1] + 13 x[7,1] + 14 x[8,1] + 20 x[9,1] + 12 x[10,1] + 20 x[11,1] + 15 x[12,1] + 11 x[13,1] + 3 x[14,1] + 8 x[15,1] + 14 x[1,2] + 14 x[2,2] + 3 x[3,2] + 8 x[4,2] + 10 x[5,2] + 11 x[6,2] + 14 x[7,2] + 18 x[8,2] + 8 x[9,2] + 6 x[10,2] + 4 x[11,2] + 14 x[12,2] + 11 x[13,2] + 12 x[14,2] + 20 x[15,2] + [[...165 terms omitted...]] + 15 x[1,14] + 16 x[2,14] + x[3,14] + 18 x[4,14] + 18 x[5,14] + 9 x[6,14] + 7 x[7,14] + 16 x[8,14] + 12 x[9,14] + 5 x[10,14] + 2 x[11,14] + 11 x[12,14] + 18 x[13,14] + 9 x[14,14] + 16 x[15,14] + 7 x[1,15] + 12 x[2,15] + 10 x[3,15] + 18 x[4,15] + 5 x[5,15] + 9 x[6,15] + 10 x[7,15] + 3 x[8,15] + 6 x[9,15] + 17 x[10,15] + 18 x[11,15] + 4 x[12,15] + 18 x[13,15] + 12 x[14,15] + 16 x[15,15] = 53.99999808424983
julia> old_rhs = normalized_rhs(c)
53.99999808424983
julia> set_normalized_rhs(c, round(Int, old_rhs))
julia> optimize!(model)
Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (mac64[arm] - Darwin 24.6.0 24G419)
CPU model: Apple M4
Thread count: 10 physical cores, 10 logical processors, using up to 10 threads
WLS license 722777 - registered to JuMP Development
Optimize a model with 33 rows, 225 columns and 1125 nonzeros (Min)
Model fingerprint: 0x82831798
Model has 225 linear objective coefficients
Variable types: 0 continuous, 225 integer (225 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+01]
Objective range [1e+01, 6e+01]
Bounds range [0e+00, 0e+00]
RHS range [1e+00, 3e+02]
Presolve removed 0 rows and 26 columns
Presolve time: 0.01s
Presolved: 33 rows, 199 columns, 914 nonzeros
Variable types: 0 continuous, 199 integer (199 binary)
Root relaxation: objective 3.344286e+02, 61 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 334.42857 0 11 - 334.42857 - - 0s
H 0 0 356.0000000 342.41176 3.82% - 0s
0 0 356.00000 0 25 356.00000 356.00000 0.00% - 0s
Cutting planes:
Gomory: 3
Cover: 3
MIR: 3
StrongCG: 1
GUB cover: 3
Zero half: 3
RLT: 2
Explored 1 nodes (66 simplex iterations) in 0.01 seconds (0.01 work units)
Thread count was 10 (of 10 available processors)
Solution count 1: 356
Optimal solution found (tolerance 1.00e-04)
Best objective 3.560000000000e+02, best bound 3.560000000000e+02, gap 0.0000%
User-callback calls 54, time in user-callback 0.00 sec
julia> set_normalized_rhs(c, old_rhs)
julia> optimize!(model)
Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (mac64[arm] - Darwin 24.6.0 24G419)
CPU model: Apple M4
Thread count: 10 physical cores, 10 logical processors, using up to 10 threads
WLS license 722777 - registered to JuMP Development
Optimize a model with 33 rows, 225 columns and 1125 nonzeros (Min)
Model fingerprint: 0x516d49ba
Model has 225 linear objective coefficients
Variable types: 0 continuous, 225 integer (225 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+01]
Objective range [1e+01, 6e+01]
Bounds range [0e+00, 0e+00]
RHS range [1e+00, 3e+02]
Loaded MIP start from previous solve with objective 356
Presolve time: 0.00s
Explored 0 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 1 (of 10 available processors)
Solution count 1: 356
Model is unbounded
Warning: max constraint violation (1.9158e-06) exceeds tolerance
Best objective 3.560000000000e+02, best bound -, gap -
User-callback calls 29, time in user-callback 0.00 sec |
|
Regarding the infeasibility, it's right on the edge so can go either way. You can generate an IIS from the model linked above: Looking at
This probably won't work in general. The previous solve found an objective value of
I can't reproduce it either - can you extract a MIP start for the second solve? Setting
|
|
Yeah no issue with the optimal/infeasible decision. Understandable and expected. The issue is the "unbounded" status. |
|
This is consistently reproducible for me: julia> using JuMP, Gurobi
julia> model = read_from_file("/tmp/model.mps")
A JuMP Model
├ solver: none
├ objective_sense: MIN_SENSE
│ └ objective_function_type: AffExpr
├ num_variables: 225
├ num_constraints: 258
│ ├ AffExpr in MOI.EqualTo{Float64}: 31
│ ├ AffExpr in MOI.LessThan{Float64}: 2
│ └ VariableRef in MOI.ZeroOne: 225
└ Names registered in the model: none
julia> set_optimizer(model, Gurobi.Optimizer)
Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID to value 722777
WLS license 722777 - registered to JuMP Development
julia> c = constraint_by_name(model, "R74")
R74 : 8 x[1,1] + 7 x[2,1] + 13 x[3,1] + 14 x[4,1] + 5 x[5,1] + 12 x[6,1] + 13 x[7,1] + 14 x[8,1] + 20 x[9,1] + 12 x[10,1] + 20 x[11,1] + 15 x[12,1] + 11 x[13,1] + 3 x[14,1] + 8 x[15,1] + 14 x[1,2] + 14 x[2,2] + 3 x[3,2] + 8 x[4,2] + 10 x[5,2] + 11 x[6,2] + 14 x[7,2] + 18 x[8,2] + 8 x[9,2] + 6 x[10,2] + 4 x[11,2] + 14 x[12,2] + 11 x[13,2] + 12 x[14,2] + 20 x[15,2] + [[...165 terms omitted...]] + 15 x[1,14] + 16 x[2,14] + x[3,14] + 18 x[4,14] + 18 x[5,14] + 9 x[6,14] + 7 x[7,14] + 16 x[8,14] + 12 x[9,14] + 5 x[10,14] + 2 x[11,14] + 11 x[12,14] + 18 x[13,14] + 9 x[14,14] + 16 x[15,14] + 7 x[1,15] + 12 x[2,15] + 10 x[3,15] + 18 x[4,15] + 5 x[5,15] + 9 x[6,15] + 10 x[7,15] + 3 x[8,15] + 6 x[9,15] + 17 x[10,15] + 18 x[11,15] + 4 x[12,15] + 18 x[13,15] + 12 x[14,15] + 16 x[15,15] = 53.99999808424983
julia> old_rhs = normalized_rhs(c)
53.99999808424983
julia> set_normalized_rhs(c, round(Int, old_rhs))
julia> optimize!(model)
Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (mac64[arm] - Darwin 24.6.0 24G419)
CPU model: Apple M4
Thread count: 10 physical cores, 10 logical processors, using up to 10 threads
WLS license 722777 - registered to JuMP Development
Optimize a model with 33 rows, 225 columns and 1125 nonzeros (Min)
Model fingerprint: 0x82831798
Model has 225 linear objective coefficients
Variable types: 0 continuous, 225 integer (225 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+01]
Objective range [1e+01, 6e+01]
Bounds range [0e+00, 0e+00]
RHS range [1e+00, 3e+02]
Presolve removed 0 rows and 26 columns
Presolve time: 0.00s
Presolved: 33 rows, 199 columns, 914 nonzeros
Variable types: 0 continuous, 199 integer (199 binary)
Root relaxation: objective 3.344286e+02, 61 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 334.42857 0 11 - 334.42857 - - 0s
H 0 0 356.0000000 342.41176 3.82% - 0s
0 0 356.00000 0 25 356.00000 356.00000 0.00% - 0s
Cutting planes:
Gomory: 3
Cover: 3
MIR: 3
StrongCG: 1
GUB cover: 3
Zero half: 3
RLT: 2
Explored 1 nodes (66 simplex iterations) in 0.01 seconds (0.01 work units)
Thread count was 10 (of 10 available processors)
Solution count 1: 356
Optimal solution found (tolerance 1.00e-04)
Best objective 3.560000000000e+02, best bound 3.560000000000e+02, gap 0.0000%
User-callback calls 50, time in user-callback 0.00 sec
julia> set_normalized_rhs(c, old_rhs)
julia> optimize!(model)
Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (mac64[arm] - Darwin 24.6.0 24G419)
CPU model: Apple M4
Thread count: 10 physical cores, 10 logical processors, using up to 10 threads
WLS license 722777 - registered to JuMP Development
Optimize a model with 33 rows, 225 columns and 1125 nonzeros (Min)
Model fingerprint: 0x516d49ba
Model has 225 linear objective coefficients
Variable types: 0 continuous, 225 integer (225 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+01]
Objective range [1e+01, 6e+01]
Bounds range [0e+00, 0e+00]
RHS range [1e+00, 3e+02]
Loaded MIP start from previous solve with objective 356
Presolve time: 0.00s
Explored 0 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 1 (of 10 available processors)
Solution count 1: 356
Model is unbounded
Warning: max constraint violation (1.9158e-06) exceeds tolerance
Best objective 3.560000000000e+02, best bound -, gap -
User-callback calls 29, time in user-callback 0.00 sec
julia> solution_summary(model)
solution_summary(; result = 1, verbose = false)
├ solver_name : Gurobi
├ Termination
│ ├ termination_status : DUAL_INFEASIBLE
│ ├ result_count : 1
│ └ raw_status : Model was proven to be unbounded. Important note: an unbounded status indicates the presence of an unbounded ray that allows the objective to improve without limit. It says nothing about whether the model has a feasible solution. If you require information on feasibility, you should set the objective to zero and reoptimize.
├ Solution (result = 1)
│ ├ primal_status : UNKNOWN_RESULT_STATUS
│ ├ dual_status : NO_SOLUTION
│ └ objective_value : 3.56000e+02
└ Work counters
├ solve_time (sec) : 1.69039e-04
├ simplex_iterations : 0
├ barrier_iterations : 0
└ node_count : 0 |
|
Thanks, I can reproduce it, but not outside of Julia/JuMP. I'll see what I can find. Any chance our MPS readers are disagreeing again (the initial model fingerprints are different)? |
Yes. JuMP annoyingly doesn't preserve row order when reading MPS files. Let me make you Gurobi's MPS: using JuMP, Gurobi
model = read_from_file("/tmp/model.mps")
set_optimizer(model, Gurobi.Optimizer)
MOI.Utilities.attach_optimizer(model)
GRBwrite(unsafe_backend(model), "/tmp/model-grb.mps") |
|
Something is very unusual. julia> model = read_from_file("/tmp/model.mps")
A JuMP Model
├ solver: none
├ objective_sense: MIN_SENSE
│ └ objective_function_type: AffExpr
├ num_variables: 225
├ num_constraints: 258
│ ├ AffExpr in MOI.EqualTo{Float64}: 31
│ ├ AffExpr in MOI.LessThan{Float64}: 2
│ └ VariableRef in MOI.ZeroOne: 225
└ Names registered in the model: none
julia> set_optimizer(model, Gurobi.Optimizer)
Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID to value 722777
WLS license 722777 - registered to JuMP Development
julia> optimize!(model)
Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (mac64[arm] - Darwin 24.6.0 24G419)
CPU model: Apple M4
Thread count: 10 physical cores, 10 logical processors, using up to 10 threads
WLS license 722777 - registered to JuMP Development
Optimize a model with 33 rows, 225 columns and 1125 nonzeros (Min)
Model fingerprint: 0x516d49ba
Model has 225 linear objective coefficients
Variable types: 0 continuous, 225 integer (225 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+01]
Objective range [1e+01, 6e+01]
Bounds range [0e+00, 0e+00]
RHS range [1e+00, 3e+02]
Presolve time: 0.00s
Explored 0 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 1 (of 10 available processors)
Solution count 0
Model is infeasible or unbounded
Best objective -, best bound -, gap -
User-callback calls 25, time in user-callback 0.00 sec
julia> GRBwrite(unsafe_backend(model), "/tmp/bug-198.mps")
0
julia> env = Gurobi.Env()
Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID to value 722777
WLS license 722777 - registered to JuMP Development
Gurobi.Env(Ptr{Nothing}(0x000000012ce4e200), false, 0)
julia> modelP = Ref{Ptr{Cvoid}}(C_NULL)
Base.RefValue{Ptr{Nothing}}(Ptr{Nothing}(0x0000000000000000))
julia> @assert GRBreadmodel(env, "/tmp/bug-198.mps", modelP) == 0
Read MPS format model from file /tmp/bug-198.mps
Reading time = 0.00 seconds
: 33 rows, 225 columns, 1125 nonzeros
julia> model = modelP[]
Ptr{Nothing}(0x000000012c3f9370)
julia> GRBoptimize(model)
Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (mac64[arm] - Darwin 24.6.0 24G419)
CPU model: Apple M4
Thread count: 10 physical cores, 10 logical processors, using up to 10 threads
WLS license 722777 - registered to JuMP Development
Optimize a model with 33 rows, 225 columns and 1125 nonzeros (Min)
Model fingerprint: 0xe6a04753
Model has 225 linear objective coefficients
Variable types: 0 continuous, 225 integer (225 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+01]
Objective range [1e+01, 6e+01]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 3e+02]
Presolve time: 0.00s
Explored 0 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 1 (of 10 available processors)
Solution count 0
Model is infeasible
Best objective -, best bound -, gap -
0
julia> GRBfree(model) |
| inner, | ||
| scalars[k], | ||
| MOI.EqualTo(zₖ), | ||
| MOI.LessThan(zₖ + 1e-5), |
There was a problem hiding this comment.
So with this and my benchmark, I now get a few extra solutions. For example, there are these mostly equivalent solutions:
julia> Y[68], Y[62]
([54.00000098502844, 128.9999926496984, 127.00000464334889], [54.0, 129.0, 127.0])I really need to go away and rethink how we handle tolerances across all of the algorithms. Things are much more complicated when we want to return multiple solutions than the hierarchical algorithm implemented by Gurobi.
|
The inconsistency is coming from Gurobi.jl defining binary variables to be "free" (this is excluded when writing an MPS file from Gurobi, it seems). This leads to wrong conclusion of unboundedness. I set bounds when adding binary variables see: jump-dev/Gurobi.jl@17aa158 this seems to work. We will also fix the inconsistency on our side. |
The bounds when writing MPS files and/or the bug in reporting unbounded? I've opened a PR to Gurobi.jl with your change. It probably needs some changes to catch the intersection with other bounds. |
|
Inspired by @torressa's PR, another way to "fix" this was to add variable bounds to the binary variables: function build_model(filename)
data = readlines(filename)
P = parse(Int, data[1])
N = parse(Int, data[2])
C = [
reduce(vcat, (parse.(Int, split(data[2+p*N+n]))' for n in 1:N))
for p in 0:(P-1)
]
model = Model()
@variable(model, 0 <= x[1:N, 1:N] <= 1, Bin)
@constraint(model, [i in 1:N], sum(x[i, :]) == 1)
@constraint(model, [j in 1:N], sum(x[:, j]) == 1)
@objective(model, Min, [sum(C[p] .* x) for p in 1:P])
return model
endThen, even without this PR, we get the expected 878 solutions: |
From a previous solve, z_k may be a value like z+eps where z in Z and eps is the feasibility tolerance. However, if scalars[k] is integer valued, then presolve may (somewhat reasonably) deduce that this problem is infeasible. Instead of rounding z_k, changing EqualTo to LessThan seemed to work on the instances I have. I don't know why. There is also a weird situation in which Gurobi declared a problem unbounded. I can't reproduce without running the entire thing, but I think it is a mix of presolve proving that there is no finite solution (because its infeasible) and yet there being a "feasible" MIP solution in memory from the previous solve. It comes down to the weird mix of tolerances in MIP starts, presolve, simplex, and this equality constraint. This is most probably a bug in Gurobi, in that it should either report optimal or infeasible. But not unbounded. I've seem something similar in SDDP.jl.
AP_p-3_n-15_ins-1.txt