Bilevel optimization

[Yark]

Hi dear professor

I tried to write a simple bilevel optimization case, but it seems that something went wrong.

When I run the program, gurobi tells me that problem = 1 and the model is infeasible.

But no matter how I adjust the parameters, it is still infeasible.

So I use the inner optimization alone as a program (like test_inner.m) and assign values to the outer variables, I get the correct result. 

What confuses me is that the value of the outer variables I gave does not exceed the [lb, ub] constraint I set. This shows that the inner model is feasible under the current outer variable values. But why is outer layer optimization always not feasible in bilevel optimization?

My code is in the attached file.

Thank you！

[Johan Löfberg]

There is no option 'usex'. You mean usex0. However, trying to initilize some values won't work here, as you haven't initialized all of them. It rarely helps the solver to only give a partial  warm-start (and on MILPs it typically doesn't speed thngs up anyway, as the hard part is not to find the optimal solution, but to prove that it is optimal)

Having said that, you call it simple_bilevel, and then you define a large bilevel program. There is no chance the naive algorithm in solvebilevel will be able to solve that. It is intended for small academic examples. The very least you have to do is to use an external algorithm 'bilevel.algorithm','external', i.e. tell solvebilevel to simply setup the KKT conditions and then throw a MILP solver on that (i.e tell it to implement the stuff done here https://yalmip.github.io/command/kkt/). Manual branching on slackness will never work on this size

>> opt = sdpsettings('bilevel.algorithm','external');

solvebilevel(CO, OO, CI, OI, y, opt)

Starting derivation of dual bounds (can be turned off using option kkt.dualbounds)

*Computing 192 primal bounds (required for dual bounds)

*Computing 288 bounds on parameterized RHS (required for dual bounds)

Parameterized RHS in equality not yet supported in lifting based bounds

*Computing 288 dual bounds

*Warning: Some of the dual upper bounds are infinite

Academic license - for non-commercial use only

Optimize a model with 1598 rows, 918 columns and 2964 nonzeros

Variable types: 624 continuous, 294 integer (288 binary)

Coefficient statistics:

  Matrix range     [1e+00, 5e+04]

  Objective range  [1e+01, 8e+04]

  Bounds range     [1e+00, 1e+00]

  RHS range        [1e-02, 5e+04]

Presolve removed 1105 rows and 504 columns

Presolve time: 0.03s

Presolved: 493 rows, 414 columns, 1520 nonzeros

Variable types: 227 continuous, 187 integer (181 binary)

Found heuristic solution: objective 2760754.9360

Found heuristic solution: objective 2513012.5380

Root relaxation: objective 2.356647e+06, 206 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work

 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 2356646.59    0    7 2513012.54 2356646.59  6.22%     -    0s

H    0     0                    2485927.4520 2356646.59  5.20%     -    0s

     0     0 2368632.53    0   22 2485927.45 2368632.53  4.72%     -    0s

     0     0 2368632.53    0   22 2485927.45 2368632.53  4.72%     -    0s

     0     0 2368632.53    0   38 2485927.45 2368632.53  4.72%     -    0s

     0     0 2368632.53    0   39 2485927.45 2368632.53  4.72%     -    0s

     0     0 2368632.53    0    3 2485927.45 2368632.53  4.72%     -    0s

H    0     0                    2371790.0620 2368632.53  0.13%     -    0s

     0     0 2370286.89    0    8 2371790.06 2370286.89  0.06%     -    0s

     0     0 2370287.19    0    4 2371790.06 2370287.19  0.06%     -    0s

     0     0 2370287.19    0    2 2371790.06 2370287.19  0.06%     -    0s

     0     0 2370305.52    0    1 2371790.06 2370305.52  0.06%     -    0s

H    0     0                    2370305.5230 2370305.52  0.00%     -    0s

Cutting planes:

  MIR: 2

Explored 1 nodes (567 simplex iterations) in 0.27 seconds

Thread count was 4 (of 4 available processors)

Solution count 5: 2.37031e+06 2.37179e+06 2.48593e+06 ... 2.76075e+06

Optimal solution found (tolerance 1.00e-04)

Best objective 2.370305523000e+06, best bound 2.370305523000e+06, gap 0.0000%

ans = 

  struct with fields:

    yalmipversion: '20190425'

       yalmiptime: 0.2534

       solvertime: 0.3596

             info: 'Successfully solved (GUROBI-GUROBI)'

          problem: 0

[Yark]

Dear Professor, thank you very much! This problem bothered me for two days. I finally know what the reason is, maybe I should be learning about kkt. Thank you again!

[Yark]

Dear Professor,

My question seems to take a long time to calculate, but I don't know how to set the timelimit and gap in the bilevel. Can you tell me how to set it up. thank!

[Johan Löfberg]

According to my display above, it is solved in a matter of seconds. You have to have a very good MILP solver installed (gurobi/cplex/mosek)

[Yark]

These codes are only part of the code. In the original program, there were many variables, and I also installed gurobi. But unfortunately, after running for 7000s, the gap is 14.7%, which is about 0.1% decrease every 5 minutes. So I want to set the time or gap to see the results.

At the same time, I found that using this bilevel optimization is not as fast as heuristic algorithm (like NSGA-II) + gurobi.

Or is there something wrong with my program?

[Johan Löfberg]

As always

https://yalmip.github.io/slowmilp

To tune the solver and various settings, you tune it via sdpsettings as always.

You should probably spend more time on creating stronger cuts for your KKT models, bounding duals etc. With your massive size, hoping a generic implementtion will work is very optimistic

[Yark]

Dear Professor

There are at most 2500 sdpvar variables in my program. After reading your reply and https://yalmip.github.io/slowmilp, I don't think it should take so long. I think you are right, so I went to see https://yalmip.github.io/command/kkt/

But I find it so simple that I don't know how to use it.

[Johan Löfberg]

The "small" problem you sent has 288 inner inequalities, so at least 288 binary variables already from the KKT representation of the inner program. After reading that link and realizing 2^288 is mind-bogglingly much larger than 2^72, you should see that it just as well can take inf years to solve it.

[Yark]

I understand. In addition, I still have a doubt, whether the bilevel optimization will ultimately get the optimal solution, or just the optimal solution for the outer objective.

Because I found that no matter how the inner layer parameters are adjusted, it does not seem to affect the final objective of the outer layer optimization.

[Johan Löfberg]

??

by defnition, bilevel optimization computes the optimal solution to the, no surprise, bilevel problem

then in practice you might of course have issues with computing this due numerical problems, computational intractability, bad bounds for duals used in big-m representation of complementarity etc etc.