`dual_objective_value` in JuMP.solution_summary

Specific Domains Optimization (Mathematical)
1

I have a question: what is the dual_objective_value below

julia> JuMP.solution_summary(model)
solution_summary(; result = 1, verbose = false)
├ solver_name          : Gurobi
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ ├ raw_status         : Model was solved to optimality (subject to tolerances), and an optimal solution is available.
│ └ objective_bound    : 3.14159e+00
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : FEASIBLE_POINT
│ ├ objective_value      : 3.14159e+00
│ └ dual_objective_value : 3.14159e+00
└ Work counters
  ├ solve_time (sec)   : 2.53916e-04
  ├ simplex_iterations : 0
  ├ barrier_iterations : 0
  └ node_count         : 0

Which Gurobi’s attribute did it refer to? (I fail to identify)

2

It doesn’t correspond to any Gurobi attribute.

dual_objective_value(model) is equivalent to querying the MOI.DualObjectiveValue attribute. This attribute returns the objective value of the dual problem.

Here’s the code:

Because Gurobi doesn’t provide access to the dual objective value, we manually compute it based on the constraint duals.

Gurobi does have ObjBound, but I think we had cases in the past where this was not always equivalent to the dual objective for LPs: Model Attributes - Gurobi Optimizer Reference Manual

1 Like
3

Inspired by a recent post where an LP is solved to suboptimality, I think it makes some sense to have relative gap/feasible solution/objective bound for LPs also. But it seems that Gurobi only has MIPGap, not applied to LPs. (no more comments unless I find that more needs are required in the future.)

I believe we can make the information of solution_summary richer.

4

I wonder what is the fastest way to generate a Benders’ feasibility cut from the 2nd-stage subproblem given an improper 1st-stage trial solution, with the Gurobi solver?

I already has a conventional method myself, which manually slack the 2nd-stage subproblem so it can be solved to OPTIMAL, and then I can generate a feasibility cut using reduced cost info and objective_value. But with my method I need to switch the status of my subproblem (so I might sometimes need to solve twice—calling optimize! twice). I wonder if retrieving information directly under INFEASIBLE status would be better. But I’m yet unclear how to do

Gurobi Optimizer version 13.0.1 build v13.0.1rc0 (linux64 - "Ubuntu 24.04.4 LTS")

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:
InfUnbdInfo  1
Threads  1

Optimize a model with 7178 rows, 11042 columns and 26012 nonzeros (Min)
Model fingerprint: 0x3bc01260
Model has 2 linear objective coefficients
Coefficient statistics:
  Matrix range     [2e-01, 2e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e-02, 1e+02]
  RHS range        [1e-01, 2e+02]

Presolve removed 7156 rows and 9054 columns
Presolve time: 0.00s
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0      handle free variables                          0s

Solved in 2278 iterations and 0.05 seconds (0.08 work units)
Infeasible model
solution_summary(; result = 1, verbose = false)
├ solver_name          : Gurobi
├ Termination
│ ├ termination_status : INFEASIBLE
│ ├ result_count       : 0
│ ├ raw_status         : Model was proven to be infeasible.
│ └ objective_bound    : 7.52770e+00
├ Solution (result = 1)
│ ├ primal_status        : NO_SOLUTION
│ └ dual_status          : INFEASIBILITY_CERTIFICATE
└ Work counters
  ├ solve_time (sec)   : 4.97351e-02
  ├ simplex_iterations : 2278
  ├ barrier_iterations : 0
  └ node_count         : 0
ERROR: Result index of attribute MathOptInterface.DualObjectiveValue(1) out of bounds. There are currently 0 solution(s) in the model.
Stacktrace:
 [1] check_result_index_bounds
   @ ~/.julia/packages/MathOptInterface/Q3V1z/src/attributes.jl:238 [inlined]
 [2] get
   @ ~/.julia/packages/Gurobi/K2XSK/src/MOI_wrapper/MOI_wrapper.jl:3376 [inlined]
 [3] _moi_get_result(model::Gurobi.Optimizer, args::MathOptInterface.DualObjectiveValue)
   @ JuMP ~/.julia/packages/JuMP/0tD10/src/optimizer_interface.jl:1191
 [4] get(model::JuMP.Model, attr::MathOptInterface.DualObjectiveValue)
   @ JuMP ~/.julia/packages/JuMP/0tD10/src/optimizer_interface.jl:1220
 [5] dual_objective_value(model::JuMP.Model; result::Int64)
   @ JuMP ~/.julia/packages/JuMP/0tD10/src/objective.jl:167
5

That doesn’t look right. Can you open an issue in Gurobi.jl with a reproducible example? The smaller the better.

6

Why?

The bigger question is

  • How to generate a feasibility cut with the Gurobi Solver, directly under INFEASIBLE status? use FarkasDual for every constraint? Should we follow this routine?

Currently I have a method myself that works good (add artificial slack variable in the 2nd-stage problem so it is always feasible). My methods has pros and cons:

  • pros: only need the primal formulation, only need to watch the linking constraint, only need the RC variable attribute to get the slope of the feasibility cut and need the ObjVal model attribute to determine the intercept of the feasibility cut.
  • cons: I need to switch the status during performing the Benders decomposition algorithm: to “turn on” and later on “turn off” the artificial variables.

I wonder what is the pros and cons of the “standard” method (as written in the JuMP’s doc, with HiGHS solver)

7

I tend to believe that results are better off retrieved directly from the solver (e.g. they use highly-optimized low-level code, so they should be the fastest).

For example, we can use JuMP.value to retrieve the numeric value of a JuMP.AffExpr, but that entails JuMP calculation, which might be slow. And It might be seen as a performance trap.

A better method should be

JuMP.@variable(m, expr)
JuMP.@constraint(m, expr == #= The long expression =#)
optimize!(m)
JuMP.value(expr) # This directly retrieve a Float64 without further calculation in julia

I think the same issue lies with JuMP.dual_objective_value, if it performs additional calculation in julia. (am I right?)

8

That doesn’t look right.

Why?

Because JuMP threw an error that there were no solutions, but there is an infeasibility certificate.

The result count should be 1 if there is a certificate:

How to generate a feasibility cut with the Gurobi Solver

If you’re doing Benders, I strongly encourage adding variables to represent the incoming state variables. It means that you need to care only about the Farkas dual of the FixRef constraints for the incoming state variables. Basically: what the tutorial does.

I think the same issue lies with JuMP.dual_objective_value, if it performs additional calculation in julia. (am I right?)

Yes, currently we haven’t implemented a method in Gurobi.jl to compute the dual objective value, so it uses a relatively expensive fallback.

9

Can we exchange our method on generating feasibility cut? Using raw info retrieved from Gurobi model/variable attributes, (I have no concrete idea how you’d like to implement dual_objective_value and reduced_cost under INFEASIBLE status). (It’s perhaps not correct to retrieve the RC variable attribute from Gurobi under a INFEASIBLE status.)

My current guess is that you need to retrieve the constraint attribute FarkasDual of all constraints (may be much more than the linking ones) in the 2nd-stage subproblem. And then calculate the intercept of the feasibility cut.

10

I’m not sure if it’s my reason: for performance reasons, I used many direct Gurobi C-API, e.g. I used Gurobi.GRBoptimize instead of JuMP.optimize!.

Now I’m used to build the physical model with JuMP, but use Gurobi C-API in the performance critical part.

11

I guess the related code are at

function _dual_objective_value(
    model::MOI.ModelLike,
    ::Type{F},
    ::Type{S},
    ::Type{T},
    result_index::Integer,
)::T where {T,F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
    value = zero(T)
    if F == variable_function_type(S) && !_variable_set_in_dual_objective(S)
        # Early return. This is a constraint like x in R_+, so no contribution
        # appears in the dual objective.
        return value
    end
    for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
        constant = MOI.constant(MOI.get(model, MOI.ConstraintFunction(), ci), T)
        set = MOI.get(model, MOI.ConstraintSet(), ci)
        dual = MOI.get(model, MOI.ConstraintDual(result_index), ci)
        value += _dual_objective_dot(constant, dual, set)
    end
    return value
end

and

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},<:Any},
)
    _throw_if_optimize_in_progress(model, attr)
    MOI.check_result_index_bounds(model, attr)
    valueP = Ref{Cdouble}()
    row = Cint(_info(model, c).row - 1)
    if model.has_infeasibility_cert
        ret = GRBgetdblattrelement(model, "FarkasDual", row, valueP)
        _check_ret(model, ret)
        return -valueP[]
    end
    ret = GRBgetdblattrelement(model, "Pi", row, valueP)
    _check_ret(model, ret)
    return _dual_multiplier(model) * valueP[]
end

and (This is in reduced_cost)

function _farkas_variable_dual(model::Optimizer, col::Cint)
    numnzP = Ref{Cint}()
    ret = GRBgetvars(model, numnzP, C_NULL, C_NULL, C_NULL, col, 1)
    _check_ret(model, ret)
    vbeg = Vector{Cint}(undef, 2)
    vind = Vector{Cint}(undef, numnzP[])
    vval = Vector{Cdouble}(undef, numnzP[])
    ret = GRBgetvars(model, numnzP, vbeg, vind, vval, col, 1)
    _check_ret(model, ret)
    λ = Vector{Cdouble}(undef, numnzP[])
    ret = GRBgetdblattrlist(model, "FarkasDual", length(vind), vind, λ)
    _check_ret(model, ret)
    return λ' * vval
end

. So I’d like to ask:

  • Do you loop over all constraints in the 2nd-stage subproblem and query the FarkasDual constraint attribute element-wise, in order to calculate the Benders’ feasibility cut?

If so, I think it’s a bit of involved.