Getting unbounded ray in JuMP

Specific Domains Optimization (Mathematical)
1

I have an unbounded linear program that I am solving using Gurobi with the attribute:

set_attribute(model,  "InfUnbdInfo", 1)

set_attribute(model, "DualReductions", 0)

and the model is indeed unbounded. I want to get an extreme ray for this unbounded model ("UnbdRay" in Gurobi), but in the JuMP manual I could not find how to do that. Any tips regarding how to get an extreme ray of an unbounded LP would be much appreciated. Also is there a solver-independent way of getting an extreme ray (provided the LP is unbounded)?

2

If primal_status(model) is INFEASIBILITY_CERTIFICATE, then value(x) is a primal unbounded ray (a certificate of dual infeasibility).

I guess I need to add an explicit example of this to the documentation.

(Edit: see [docs] add section on unbounded rays by odow · Pull Request #3945 · jump-dev/JuMP.jl · GitHub)

3

Okay that works, thanks @odow!

1 Like
4

In following up on this thread and in connection with Stuck with infeasibility certificates in JuMP.jl, I would like to understand why different solvers produce inconsistent values of ObjectiveValue for the following code:

using JuMP, Gurobi

# Data
e = 3.0
B = [[6.0, 2.0], [6.0, 2.0]]
N_B = [length(B[i]) for i in 1:2]
ȳ = [5.0, 1.0]
c = 3.0

z̄ = zeros(length(N_B), maximum(N_B)) 
π̄ = zeros(length(N_B), maximum(N_B))
σ̄ = zeros(length(N_B), maximum(N_B))
π̄ⁱⁿᶠ = zeros(length(N_B), maximum(N_B))
σ̄ⁱⁿᶠ = zeros(length(N_B), maximum(N_B))

for i in 1:length(N_B)
    model = Model(Gurobi.Optimizer)
    set_optimizer_attribute(model, "InfUnbdInfo", 1)
    @variable(model, z[i, 1:N_B[i]] >= 0)
    @objective(model, Min, sum(e * z[i, j] for j in 1:N_B[i]))
    Const = @constraint(model, [j in 1:N_B[i]], ȳ[i] + z[i, j] >= B[i][j])
    ConstBound = @constraint(model, [j in 1:N_B[i]], c >= z[i, j])
    optimize!(model)
    if termination_status(model) == MOI.OPTIMAL
        z̄[i, 1:N_B[i]] = vec(value.(z[i, :]))
        π̄[i, 1:N_B[i]] = [dual(Const[j]) for j in 1:N_B[i]]
        σ̄[i, 1:N_B[i]] = [dual(ConstBound[j]) for j in 1:N_B[i]]
    elseif termination_status(model) == MOI.INFEASIBLE
        z̄[i, 1:N_B[i]] = vec(value.(z[i, :]))
        π̄ⁱⁿᶠ[i, 1:N_B[i]] = dual.(Const)
        σ̄ⁱⁿᶠ[i, 1:N_B[i]] = [dual(ConstBound[j]) for j in 1:N_B[i]]
    end
end

The infeasibility occurs for i = 2 and j = 1, and I got the following values for z̄ when using different solvers:
Gurobi:
[1.0 0.0; 5.0 1.0],
GLPK:
[1.0 0.0; 5.0 1.0],
HiGHS:
[1.0 0.0; 1.13468e-311 2.0e-323],
and Mosek:
[1.0 0.0; 0.0 0.0].
It seems like different solvers handle infeasibility differently, leading to inconsistent values. My question is: Do I need to explicitly write the dual problem and use Farkas’ Lemma to get the dual rays for infeasible primal problems, independent of the solvers, to make the code more general and consistent across different solvers?

I’d really appreciate it if someone could explain this in the context of this. Thanks!

5

Can anyone help with this?

6

In the if termination_status(model) == MOI.INFEASIBLE branch you are querying primal value vec(value.(z[i, :])) when the primal_status(model) == NO_SOLUTION .

The result is undefined when there is an infeasibility certificate, so the difference between solvers is expected behavior. See: Infeasibility certificates · JuMP

If you are interested in feasibility cuts, see Benders decomposition · JuMP