# Getting unbounded ray in JuMP

**URL:** <https://discourse.julialang.org/t/getting-unbounded-ray-in-jump/126097>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [February 20, 2025, 2:43am UTC](https://discourse.julialang.org/t/getting-unbounded-ray-in-jump/126097 "2025-02-20T02:43:30Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Shuvomoy\_Das\_Gupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shuvomoy_das_gupta/32/10069_2.png) [@Shuvomoy\_Das\_Gupta](https://discourse.julialang.org/u/Shuvomoy_Das_Gupta)\
**Post date:** [February 20, 2025, 2:43am UTC](https://discourse.julialang.org/t/getting-unbounded-ray-in-jump/126097/1 "2025-02-20T02:43:31Z")

</div>

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

```julia
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)?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [February 20, 2025, 3:57am UTC](https://discourse.julialang.org/t/getting-unbounded-ray-in-jump/126097/2 "2025-02-20T03:57:02Z")

</div>

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](https://github.com/jump-dev/JuMP.jl/pull/3945))

---

<div class="post-metadata">

**Author:** ![Shuvomoy\_Das\_Gupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shuvomoy_das_gupta/32/10069_2.png) [@Shuvomoy\_Das\_Gupta](https://discourse.julialang.org/u/Shuvomoy_Das_Gupta)\
**Post date:** [February 20, 2025, 4:20am UTC](https://discourse.julialang.org/t/getting-unbounded-ray-in-jump/126097/3 "2025-02-20T04:20:22Z")

</div>

Okay that works, thanks @odow!

---

<div class="post-metadata">

**Author:** ![ashefa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ashefa/32/2525_2.png) [@ashefa](https://discourse.julialang.org/u/ashefa)\
**Post date:** [February 26, 2025, 11:41am UTC](https://discourse.julialang.org/t/getting-unbounded-ray-in-jump/126097/4 "2025-02-26T11:41:20Z")

</div>

In following up on this thread and in connection with [Stuck with infeasibility certificates in JuMP.jl](https://discourse.julialang.org/t/stuck-with-infeasibility-certificates-in-jump-jl/102210), I would like to understand why different solvers produce inconsistent values of [`ObjectiveValue`](https://jump.dev/MathOptInterface.jl/stable/reference/models/#MathOptInterface.ObjectiveValue) for the following code:

```julia
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](https://jump.dev/JuMP.jl/dev/moi/manual/solutions/#Primal-dual-convex-solver). Thanks!

---

<div class="post-metadata">

**Author:** ![ashefa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ashefa/32/2525_2.png) [@ashefa](https://discourse.julialang.org/u/ashefa)\
**Post date:** [February 28, 2025, 8:31am UTC](https://discourse.julialang.org/t/getting-unbounded-ray-in-jump/126097/5 "2025-02-28T08:31:07Z")

</div>

Can anyone help with this?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [February 28, 2025, 6:33pm UTC](https://discourse.julialang.org/t/getting-unbounded-ray-in-jump/126097/6 "2025-02-28T18:33:29Z")

</div>

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](https://jump.dev/JuMP.jl/dev/moi/background/infeasibility_certificates/#Infeasible-problems)

If you are interested in feasibility cuts, see [Benders decomposition · JuMP](https://jump.dev/JuMP.jl/dev/tutorials/algorithms/benders_decomposition/#Feasibility-cuts)