Dear developer,
In the past period I have significant problems in finding a result for my application and developed the following short piece of code that reproduce my issue.
The problem described above represent the case of two people needing fruits (i.e. 3 fruit the first and 6 the second one) that cost 10 cents to the global market. However, they could also buy and sell fruits between them at price p.
m_P/N represent the fruits exchanged between the users and p_N are the fruits bought from the external market by each user.
Since no user produce fruits, the only solution of this problem is that each user buys the requested fruits from the market.
Therefore: the only possible solution is m_P/N = 0, p_N = [3, 6]. That holds at any price p since users do not produce fruits.
The problem can be represented as
(Upper)
Max sum(0.01*m_N[i] + 0.003*m_P[i] for i in 1:2)
s.t.
sum(m_P[i] - m_N[i] for i in 1:2) == 0
0 <= p <= 0.1
(Lower)
s.t
max sum( - p_N[i]*0.1 - m_N[i]*(p + 0.01) + m_P[i]*(p - 0.003) for i in 1:2)
s.t.
-p_N[i] + m_P[i] - m_N[i] == - 3*i for i=1:2
0 <= p/m_N/P <= 10
The proposed code is aimed to represent the problembelow works and find the solution m_P/N = 0, p_N = [3, 6] correctly. However, the solution comes with the value of p = about 0.09.
Theoretically, this problem should be working for every value of p, as explained above, however, when I change the bounds of p, infeasibilities occur.
using JuMP, BilevelJuMP, Gurobi
m = BilevelModel(Gurobi.Optimizer, mode = BilevelJuMP.FortunyAmatMcCarlMode(primal_big_M = 1000, dual_big_M = 1000))
@variable(Upper(m), 0 <= p <= 0.1)
@variable(Lower(m), 0 <= p_N[1:2] <= 10)
@variable(Lower(m), 0 <= m_P[1:2] <= 10)
@variable(Lower(m), 0 <= m_N[1:2] <= 10)
@constraint(Upper(m), con_a, sum(m_P[i] - m_N[i] for i in 1:2) == 0)
@objective(Upper(m), Max, sum(0.01*m_N[i] + 0.003*m_P[i] for i in 1:2))
@constraint(Lower(m), con_b[i=1:2], - p_N[i] + m_P[i] - m_N[i] == - 3*i)
@objective(Lower(m), Max, sum( - p_N[i]*0.1 - m_N[i]*(p + 0.01) + m_P[i]*(p - 0.003) for i in 1:2))
optimize!(m, bilevel_prob = "prova.lp", problem_format = MOI.FileFormats.FORMAT_AUTOMATIC)
_p = value.(p)
_p_N = value.(p_N)
_m_P = value.(m_P)
_m_N = value.(m_N)
@show _p
@show _p_N
@show _m_P
@show _m_N
For example, if after the previous code, these lines are added:
set_upper_bound(p, 0.05)
optimize!(m, bilevel_prob = "prova2.lp") # error
The new problem becomes infeasible, which is very weird.
To confirm that the situation is weird, if also the bounds of variables m_P/N are set to 0, then the optimization works.
However, as stated in the beginning, m_P/N==0 was a result of the first optimization, so the fact that by imposing these additional constraints, the problem works suggest that there might be a bug.
for i = 1:2
set_upper_bound(m_P[i], 0)
set_upper_bound(m_N[i], 0)
end
set_upper_bound(p, 0.05)
optimize!(m, bilevel_prob = "prova3.lp") # works but the solution hasn't changed
I am trying to investigate the issue, because I believe that this hidden bug may have significant silent effects in the optimality of the proposed problems.
P.S. I found it more useful to have the lp format. To do so, I changed the function print_lp as follows
function print_lp(m, name, problem_format = MOI.FileFormats.FORMAT_AUTOMATIC)
dest = MOI.FileFormats.Model(format = problem_format, filename = name)
MOI.copy_to(dest, m)
MOI.write_to_file(dest, name)
end
File containing all code
test_bug.txt