The output of Gurobi's solution isn't qualified as its own input in terms of feasibility

I write decomposition algorithms with the Gurobi solver, so occasionally I need to use Gurobi’s solution as data input to another model, which is also solved by Gurobi. But I find it annoying that the “$THE_TITLE”.

julia> JuMP.owner_model(X[2583].bU_1[29])
A JuMP Model
├ mode: DIRECT
├ solver: Gurobi
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 1080
├ num_constraints: 2025
│ ├ JuMP.AffExpr in MOI.GreaterThan{Float64}: 33
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 384
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 552
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 360
│ ├ JuMP.VariableRef in MOI.ZeroOne: 336
│ └ JuMP.VariableRef in MOI.LessThan{Float64}: 360
└ Names registered in the model
  └ :primobj

julia> JuMP.owner_model(X[2583].bU_1[29]) === inn[2583]
true

julia> let v = X[2583].bU_1[29], m = inn[2583]
           @assert JuMP.is_binary(v)
           JuMP.optimize!(m)
           @assert JuMP.termination_status(m) == JuMP.OPTIMAL
           @assert JuMP.primal_status(m) == JuMP.FEASIBLE_POINT
           println(JuMP.value(v))
       end
Gurobi Optimizer version 12.0.3 build v12.0.3rc0 (linux64 - "Debian GNU/Linux 12 (bookworm)")

CPU model: AMD EPYC 7763 64-Core Processor, instruction set [SSE2|AVX|AVX2]
Thread count: 128 physical cores, 256 logical processors, using up to 1 threads

Non-default parameters:
TimeLimit  300
Threads  1

Optimize a model with 777 rows, 1080 columns and 3959 nonzeros
Model fingerprint: 0x43e59afa
Variable types: 744 continuous, 336 integer (336 binary)
Coefficient statistics:
  Matrix range     [2e-01, 1e+01]
  Objective range  [3e+00, 8e+00]
  Bounds range     [8e+00, 5e+01]
  RHS range        [1e+00, 1e+02]
Presolved: 595 rows, 802 columns, 2940 nonzeros

Continuing optimization...


Cutting planes:
  Gomory: 3
  Cover: 1
  MIR: 30
  Flow cover: 7

Explored 16 nodes (1484 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 1 (of 256 available processors)

Solution count 10: 10214.6 10215.7 10215.9 ... 11370.6

Optimal solution found (tolerance 1.00e-04)
Best objective 1.021463213362e+04, best bound 1.021463213362e+04, gap 0.0000%

User-callback calls 37, time in user-callback 0.00 sec
3.1340395433073535e-6

The annoying consequence is that

julia> m = JuMP.Model(Gurobi.Optimizer)

julia> JuMP.@variable(m, x, Int)

julia> JuMP.@constraint(m, x == 3.134e-6)

julia> JuMP.optimize!(m)
Gurobi Optimizer version 12.0.3 build v12.0.3rc0 (linux64 - "Debian GNU/Linux 12 (bookworm)")

Model is infeasible
Best objective -, best bound -, gap -

I guess we need to “rectify” the raw solution from Gurobi somehow…

Hi Walter,
Is there a question there somewhere?

To let us know the issue so we can develop better software.

By the way, one thing I love about JuMP is that it has very high-quality docs, particularly comparing to the julia’s. I agree that julia the software is of very high quality, but I think it’s doc has room for improvement.

To let us know the issue

There’s no issue. This is expected behaviour.

The Optimal solution found (tolerance 1.00e-04) refers to the default MIP gap (Parameter Reference - Gurobi Optimizer Reference Manual), not the integer feasibility tolerance.

I think you overlooked my this line, which indicates that Gurobi deem it feasible (implies integrality feasibility.) which contradicts itself (the latter model reports INFEASIBLE.)


To a larger extent,

Even if (we assume temporarily) there’s no issue with Gurobi, I think there is room for improvement on how to overcome this difficulty with the julia language, e.g. multiple dispatch, etc.

As I said, I’m not expecting anything. (I’m not demanding anyone to reply me fast, or reply at all)
I post because of there’s room for improvement.

This is expected behaviour. Read this part of the tutorial: Tolerances and numerical issues · JuMP

Note that your second model is not exactly the same because it fixes the value of x with an affine constraint. Here’s a model that’s closer:

julia> using JuMP, Gurobi

julia> begin
           model = Model(Gurobi.Optimizer)
           @variable(model, x == 3.134e-6, Bin)
           optimize!(model)
           solution_summary(model)
       end
Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID to value 722777
WLS license 722777 - registered to JuMP Development
Gurobi Optimizer version 12.0.3 build v12.0.3rc0 (mac64[arm] - Darwin 24.6.0 24G90)

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 0 rows, 1 columns and 0 nonzeros
Model fingerprint: 0x5e1b6d3d
Variable types: 0 continuous, 1 integer (1 binary)
Coefficient statistics:
  Matrix range     [0e+00, 0e+00]
  Objective range  [0e+00, 0e+00]
  Bounds range     [3e-06, 3e-06]
  RHS range        [0e+00, 0e+00]
Presolve removed 0 rows and 1 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 1 (of 10 available processors)

Solution count 2: 0 0 

Optimal solution found (tolerance 1.00e-04)
Warning: max bound violation (3.1340e-06) exceeds tolerance
Best objective 0.000000000000e+00, best bound 0.000000000000e+00, gap 0.0000%

User-callback calls 91, time in user-callback 0.00 sec
solution_summary(; result = 1, verbose = false)
├ solver_name          : Gurobi
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 2
│ ├ raw_status         : Model was solved to optimality (subject to tolerances), and an optimal solution is available.
│ └ objective_bound    : 0.00000e+00
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ ├ objective_value      : 0.00000e+00
│ └ relative_gap         : 0.00000e+00
└ Work counters
  ├ solve_time (sec)   : 1.03211e-03
  ├ simplex_iterations : 0
  ├ barrier_iterations : 0
  └ node_count         : 0

Regardless: there are models for which two solvers (or even two algorithms within a solver) can disagree about whether the solution is feasible.

1 Like

Using fixing constraint makes some sense, but it’s not the case in my current situation so I still need to “round” the solution.

The example there is

M = 1e6
m = Model()
@variable(m, x >= 0)
@variable(m, y, Bin)
@constraint(m, x <= (M)y)

The solution could be (1, 1e-6).
In practice this solution can be interpreted as “Having a maximum yield of 1000000 if the machine is on”. So rounding the solution to (1, 0) means that we still has a yield of 1 when the machine is at off state. I think this rounding action is forgivable because 1 is really small compared to the full yield 1e6.