# 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:** 1
**Showing post:** 4

<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!

---

_[View the full topic](https://discourse.julialang.org/t/getting-unbounded-ray-in-jump/126097)._