Weird bug, maybe problem somewhere in build_bilevel

[davide-f]

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

[joaquimg]

Feel free to open a PR to improve printing.

[joaquimg]

Did you check is the formulation is correct with the LP printing?

Did you try it for a very small instance? Those are usually easier to debug.

[davide-f]

The attached one is the smallest problem that reproduces the problem and makes sense. Most of the listing makes sense, but I had to check it deeper and also the dualization process, to check whether there might be a sign or something in the wrong place. Unlikely I don't have a solution yet, I posted as soon as I had a short code able to reproduce the problem so to start a preliminary discussion so that maybe with 2N eyes we can find the solution faster.

[joaquimg]

Sure thing, I will also look at the example.

[davide-f]

In order to argument the finding, I am providing the following attached files: test_bug is the code (basically the same as the previous comment), prova_lp.txt is the "prova1.lp" file calculated by the test_bug code, and bug_check.pdf is the document that describes the calculation of the dual for the proposed example and details where the bug is.

If the calculations are correct, the bug should be in the sign of the c vector used in the dual reformulation and the inequality constraint. Probably the reference Max/Min may be reversed.

test_bug.txt prova_lp.txt test_bug.txt

[davide-f]

After I mentioned the comment in dualization package, they clarified that the convention used is different from the textbook one, and the dual reformulation is consistent with the goal; however that doesn't solve this bug... unlikely... Any ideas?

[joaquimg]

Hi @davide-f , I am taking a deeper look into this after a busy month. My previous comment was wrong, but I am taking a look anyway.

[joaquimg]

Ok, got it. The issue is a bit subtle. Everything is related to the constraint of the UPPER level:

@constraint(Upper(m), con_a, sum(m_P[i] - m_N[i] for i in 1:2) == 0)

I will take an extreme example. Take p = 0.0 and let's look at the solution second level problem.

max sum( -p_N[i]*0.1 -m_N[i]*0.01 + m_P[i]*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[i] <= 10

The actual optimal solution for this problem is:

m_P[i] = 0
m_N[1] = 3
m_N[2] = 6
p_N[i] = 0

Because m_N is cheaper than p_N. This is the only KKT point for p = 0 all other solutions are invalid.

However, this optimal solution is invalidated by the constraint:

sum(m_P[i] - m_N[i] for i in 1:2) == 0

Hence, the upper level has to find a different value of p so that the lower level can be solved to optimality.

Consider now p = 0.089

max sum( - p_N[i]*0.1 - m_N[i]*0.099 + m_P[i]*0.086 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[i] <= 10

Again, the optimal solution for this problem is:

m_P[i] = 0
m_N[1] = 3
m_N[2] = 6
p_N[i] = 0

Because m_N is cheaper than p_N. This is the only KKT point for p = 0 all other solutions are invalid.

Hence, there is no p in [0, 0.09) for which the optimal solutions of the second level are accepted. Therefore, the bilevel program is indeed infeasible whenever p is constrained to be in a closed subset of [0, 0.09).

[davide-f]

Thank you very much!