Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Revise status logic #519

Closed
torressa opened this issue Oct 2, 2023 · 9 comments
Closed

Revise status logic #519

torressa opened this issue Oct 2, 2023 · 9 comments

Comments

@torressa
Copy link
Collaborator

torressa commented Oct 2, 2023

From discourse discussion

The attribute SolCount is used in a few places in the MOI wrapper to set the status (e.g. src/MOI_wrapper/MOI_wrapper.jl#L2849).
When terminating early this may be 0 and no solution will appear available when there may be one.

@odow
Copy link
Member

odow commented Oct 2, 2023

I'm not sure I understand the problem. From the example in discourse:

julia> using JuMP

julia> using Gurobi

julia> using JSON

julia> data = JSON.parse("""
       {
           "plants": {
               "Seattle": {"capacity": 350},
               "San-Diego": {"capacity": 600}
           },
           "markets": {
               "New-York": {"demand": 300},
               "Chicago": {"demand": 300},
               "Topeka": {"demand": 300}
           },
           "distances": {
               "Seattle => New-York": 2.5,
               "Seattle => Chicago": 1.7,
               "Seattle => Topeka": 1.8,
               "San-Diego => New-York": 2.5,
               "San-Diego => Chicago": 1.8,
               "San-Diego => Topeka": 1.4
           }
       }
       """)
Dict{String, Any} with 3 entries:
  "plants"    => Dict{String, Any}("Seattle"=>Dict{String, Any}("capacity"=>350), "San-Diego"=>Dict{String, Any}("capacity"=>600))
  "distances" => Dict{String, Any}("San-Diego => New-York"=>2.5, "Seattle => Topeka"=>1.8, "Seattle => New-York"=>2.5, "San-Diego =>…
  "markets"   => Dict{String, Any}("Chicago"=>Dict{String, Any}("demand"=>300), "Topeka"=>Dict{String, Any}("demand"=>300), "New-Yor

julia> P = keys(data["plants"])
KeySet for a Dict{String, Any} with 2 entries. Keys:
  "Seattle"
  "San-Diego"

julia> M = keys(data["markets"])
KeySet for a Dict{String, Any} with 3 entries. Keys:
  "Chicago"
  "Topeka"
  "New-York"

julia> distance(p::String, m::String) = data["distances"]["$(p) => $(m)"]
distance (generic function with 1 method)

julia> model = direct_model(Gurobi.Optimizer())
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: DIRECT
Solver name: Gurobi

julia> @variable(model, x[P, M] >= 0)
2-dimensional DenseAxisArray{VariableRef,2,...} with index sets:
    Dimension 1, ["Seattle", "San-Diego"]
    Dimension 2, ["Chicago", "Topeka", "New-York"]
And data, a 2×3 Matrix{VariableRef}:
 x[Seattle,Chicago]    x[Seattle,Topeka]    x[Seattle,New-York]
 x[San-Diego,Chicago]  x[San-Diego,Topeka]  x[San-Diego,New-York]

julia> @constraint(model, [p in P], sum(x[p, :]) <= data["plants"][p]["capacity"])
1-dimensional DenseAxisArray{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape},1,...} with index sets:
    Dimension 1, ["Seattle", "San-Diego"]
And data, a 2-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape}}:
 x[Seattle,Chicago] + x[Seattle,Topeka] + x[Seattle,New-York]  350
 x[San-Diego,Chicago] + x[San-Diego,Topeka] + x[San-Diego,New-York]  600

julia> @constraint(model, [m in M], sum(x[:, m]) >= data["markets"][m]["demand"])
1-dimensional DenseAxisArray{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.GreaterThan{Float64}}, ScalarShape},1,...} with index sets:
    Dimension 1, ["Chicago", "Topeka", "New-York"]
And data, a 3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.GreaterThan{Float64}}, ScalarShape}}:
 x[Seattle,Chicago] + x[San-Diego,Chicago]  300
 x[Seattle,Topeka] + x[San-Diego,Topeka]  300
 x[Seattle,New-York] + x[San-Diego,New-York]  300

julia> @objective(model, Min, sum(distance(p, m) * x[p, m] for p in P, m in M));

julia> set_optimizer_attribute(model,"Method",2)
Set parameter Method to value 2

julia> set_optimizer_attribute(model,"Crossover",0)
Set parameter Crossover to value 0

julia> thr = 1e3
1000.0

julia> c = Float64[]
Float64[]

julia> function absolute_stopping_criterion(cb_data,cb_where)
           objbound =Ref{Cdouble}()
           Gurobi.GRBcbget(cb_data,cb_where,Gurobi.GRB_CB_BARRIER_DUALOBJ,objbound)
           objbound = objbound[]
           push!(c,objbound)
           if maximum(c) > thr
               GRBterminate(backend(model).inner)
           end
           return
       end
absolute_stopping_criterion (generic function with 1 method)

julia> MOI.set(model,Gurobi.CallbackFunction(),absolute_stopping_criterion)

julia> optimize!(model)
Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Core(TM) i5-8259U CPU @ 2.30GHz
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 5 rows, 6 columns and 12 nonzeros
Model fingerprint: 0x7b3178d3
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 2e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [3e+02, 6e+02]
Presolve time: 0.00s
Presolved: 5 rows, 6 columns, 12 nonzeros
Ordering time: 0.00s

Barrier statistics:
 AA' NZ     : 6.000e+00
 Factor NZ  : 1.500e+01
 Factor Ops : 5.500e+01 (less than 1 second per iteration)
 Threads    : 1

                  Objective                Residual
Iter       Primal          Dual         Primal    Dual     Compl     Time
   0   5.10364330e+03  0.00000000e+00  8.67e+02 0.00e+00  7.53e+02     0s
   1   2.01271371e+03  6.88171313e+02  7.54e+01 0.00e+00  1.30e+02     0s
   2   1.72708744e+03  1.63134933e+03  0.00e+00 0.00e+00  8.70e+00     0s

Barrier performed 2 iterations in 0.00 seconds (0.00 work units)
Optimization interrupted


User-callback calls 52, time in user-callback 0.02 sec

julia> valueP = Ref{Cint}()
Base.RefValue{Int32}(1195440432)

julia> ret = GRBgetintattr(unsafe_backend(model), "SolCount", valueP)
0

julia> valueP
Base.RefValue{Int32}(0)

julia> valueP[]
0

Gurobi did not find a solution, so SolCount is reported as 0.

Do you want Gurobi to report that SolCount is 1?

@torressa
Copy link
Collaborator Author

torressa commented Oct 2, 2023

SolCount is 0 as the solver has been interrupted, but a feasible solution is available (as it can be queried with the "X" attribute). I guess removing the use of SolCount to return solution values would be my preferred approach.

Even though SolCount is 0 you can access the solution via

valuePdouble = Ref{Cdouble}()
ret = GRBgetdblattrelement(unsafe_backend(model), "X", 0, valuePdouble)
println("X = $(valuePdouble[]) ret = $(ret)")

@odow
Copy link
Member

odow commented Oct 3, 2023

Even though SolCount is 0 you can access the solution via

@simonbowly is this a bug that should be fixed internally?

If Gurobi tells us that SolCount is 0, I don't think we should support querying the solution. Even though GRBgetdblattrelement returns something for this example, I don't know that it holds in general.

Have you verified that the solution is a feasible primal solution?

@simonbowly
Copy link
Collaborator

We're discussing if we can improve this. One point that came up is that intermediate LP solutions are quite rarely primal feasible. Also, a user can query the X attribute (or BarX) regardless of feasibility if a barrier solve is terminated early. On it's own this doesn't say anything about feasibility, so you should also check the solution quality attributes to verify if the point is feasible.

So I don't think it's a Gurobi bug or Gurobi.jl bug; the feasibility check is just conservative. I wouldn't suggest trying anything tricky like mapping a tolerance check to MOI.FEASIBLE_POINT in such cases. Direct C API calls to extract the relevant attributes are the way to go if really needed, but the user needs to be aware that they may not be querying a feasible point.

@torressa
Copy link
Collaborator Author

torressa commented Oct 3, 2023

One point that came up is that intermediate LP solutions are quite rarely primal feasible.

That's a very good point.

Cool thanks both. Users can query the solution via the C calls and check the quality attribute(s) (e.g. MaxVio or ComplVio, for more see here).

For future reference, the solution and quality attributes can be queried as:

# Extract variable `var` values via querying "X" attribute
for i in I
    valueP = Ref{Cdouble}()
    col = Cint(backend(model).variable_info[var[i].index].column - 1)
    ret = GRBgetdblattrelement(unsafe_backend(model), "X", col, valueP)
    println("variable = $(var[i]) has value $(valueP[])")
end

# Extract and check quality attributes
maxvio = MOI.get(model, Gurobi.ModelAttribute("MaxVio"))
complvio = MOI.get(model, Gurobi.ModelAttribute("ComplVio"))
println("maxvio=$(maxvio) complvio=$(complvio)")

@torressa torressa closed this as completed Oct 3, 2023
@odow
Copy link
Member

odow commented Oct 4, 2023

@simonbowly is there a reason that https://www.gurobi.com/documentation/current/refman/bestobjstop.html#parameter:BestObjStop applies only to MIP models?

@simonbowly
Copy link
Collaborator

The solution strategy is quite different for continuous vs MIP models. For MIPs, the solver almost always has a feasible incumbent during the search, and the user can control the search emphasis with parameters like MIPFocus and Heuristics. So in terms of termination criteria, MIPGap and BestObjStop make sense in the context of MIP.

For continuous models, the solver doesn't store feasible primal solutions it finds along its search path (if any). It's typical to have some primal and some dual infeasibility right up until the model is solved. So termination criteria which depend on having a feasible primal solution would trigger very rarely. Maybe this is a historical thing, but the approach to LP in particular has always been optimal-or-bust.

@odow
Copy link
Member

odow commented Oct 5, 2023

So termination criteria which depend on having a feasible primal solution would trigger very rarely. Maybe this is a historical thing, but the approach to LP in particular has always been optimal-or-bust.

I've had a surprising number of people ask me about this recently. It came up in discussions with the HiGHS folks as well.

It seems like if at any point you have a feasible primal solution during the solve then you could return if the objective is better than BestObjStop. You don't necessarily need to return a feasible dual solution as well.

It doesn't matter if it triggers very rarely. It seems a lot of people are starting to solve big LP power system models, and they're often fine with sub-optimal solutions, provided that they are feasible.

@simonbowly
Copy link
Collaborator

I've had a surprising number of people ask me about this recently. It came up in discussions with the HiGHS folks as well.

Thanks, interesting to know!

It seems like if at any point you have a feasible primal solution during the solve then you could return if the objective is better than BestObjStop. You don't necessarily need to return a feasible dual solution as well.
It doesn't matter if it triggers very rarely. It seems a lot of people are starting to solve big LP power system models, and they're often fine with sub-optimal solutions, provided that they are feasible.

Agreed. I understood your question as "why is it currently this way?", and I think the reason is that LP and MIP historically used different strategies and users had different requirements for the two problem types. It seems that the user requirements are converging a bit, for sure.

In the meantime, as discussed on this thread, it's currently not too hard to achieve this via a callback.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Development

No branches or pull requests

3 participants