An LP Problem using JuMP and gurobi in Julia

Question

using JuMP
using Gurobi
using MAT
using LinearAlgebra

#p_max,p_min, z, W ,b ,D are given matrix and can guarantee the correctness 

#introduce u W D b
file_u = matopen(raw"C:\Users\minboli\Desktop\matlab data\Xincode\Xincode\u.mat")
z = read(file_u,"u")

file_W = matopen(raw"C:\Users\minboli\Desktop\matlab data\Xincode\Xincode\W.mat")
W = read(file_W,"W")

file_b = matopen(raw"C:\Users\minboli\Desktop\matlab data\Xincode\Xincode\b.mat")
b = read(file_b,"b")

file_D = matopen(raw"C:\Users\minboli\Desktop\matlab data\Xincode\Xincode\sparse_matrix_D.mat")
D = read(file_D,"D")

n = length(z)
m = length(b)

function sub_problem(p_max, p_min, n, m, z, b, W, D) 
    model = Model(Gurobi.Optimizer)

    @variable(model, lambda[1:n] >= 0)
    @variable(model, mu[1:m])
    @variable(model, delta[1:m], Bin)  
    @variable(model, p[1:m])

    M = 1e5 

    @objective(model, Min, dot(lambda, z) + dot(mu, p - b))

    @constraint(model, W' * lambda + D' * mu .== 0)

  
    @constraint(model, p .<= p_max .- (1 .- delta) * M)
    @constraint(model, p .>= p_min .+ delta * M)
    @constraint(model, mu .<= M * delta)
    @constraint(model, mu .>= -M * (1 .- delta))

    optimize!(model)
    status = termination_status(model)
    if status != MOI.OPTIMAL
        error("The model did not solve correctly, status: $(status)")
    end

    return value.(lambda), value.(mu)
end

p_max_revised = reshape(p_max,:,1)
p_min_revised = reshape(p_min,:,1)

print(sub_problem(p_max_revised, p_min_revised, n, m, z, b, W, D))

After I do this optimization problem, Gurobi returns

Model is infeasible or unbounded
Best objective -, best bound -, gap -

User-callback calls 43, time in user-callback 0.00 sec
ERROR: The model did not solve correctly, status: INFEASIBLE_OR_UNBOUNDED

Infeasible or unbounded is acceptable since my optimization have encountered this kind of situation, I wonder if optimization is failed, can Gurobi return variables that do not satisfy the constraints?

I have tried to check my code, it can fulfill the objective function and subject.

Answer 1

JuMP has a number of tools to debug infeasible and unbounded models.

See the debugging tutorial:

https://jump.dev/JuMP.jl/stable/tutorials/getting_started/debugging/#Debugging-an-infeasible-model

As well as JuMP.relax_with_penalty!:

https://jump.dev/JuMP.jl/stable/api/JuMP/#JuMP.relax_with_penalty!

And tools to compute the IIS:

https://jump.dev/JuMP.jl/stable/manual/solutions/#Conflicts

Answer 2

I will address here the given question "can Gurobi return variables that do not satisfy the constraints" which is I believe your main point

Consider the following optimization problem:

x - 2y -> Max
3x >= 9
2y >= 8
x + y <= 6
x, y >= 0

Which variable in your opinion are violating the constraints?

You can immediately see that it usually makes more sense to look at constraints rather than variables. However, obviously constraints are contradictory as a set - not as a single one.

Now regardless whether you want to look at constraints or variables, you can find out "which one is violated" by introducing an artificial variable to the model. It can not be done automatically as you as the model author need to decide what type violation you want to consider.

For an example for the model above you could do:

@variable(model, x >= 0)
@variable(model, y >= 0)
@variable(model, s[1:3] >= 0)
@constraint(model, 3x + s[1] >= 9)
@constraint(model, 2y + s[2] >= 8)
@constraint(model, x + y - s[3] <= 6)

Now you have an artificial variable that describes how much model is violated (here it is continues, sometimes you might want to have a binary one multiplied by M - this is one of many decisions to be made and hence why the process cannot be automatic).

All you need to do now is to include this artificial variable in the objective:

@objective(model, Max, x - 2y - sum(M .* s))

Looking at non-zero s values at the optimal solution you will find out which constraints need to be lifted in order for the problem to be feasible.