Bilevel programming
YALMIP has built-in support for definition, setup, and solution of bilevel programming problems. The code here concentrates on the built-in solver for bilevel problems. You can of course set them up yourself, by manually deriving the KKT conditions and solving them using various techniques in YALMIP, or by using YALMIPs high-level kkt operator, as illustrated in the bilevel example.
For an introduction to bilevel optimization, see Bard 1999.
KKT conditions in bilevel programming
The class of bilevel problems that can be adressed natively by YALMIP has to have the following leader-follower (outer-inner) structure
\[\begin{aligned} \text{minimize } & f(x,y^{\star})\\ \text{subject to } & (x,y^{\star}) \in \mathcal{C}\\ &\begin{aligned} y^{\star} = & \arg \min \frac{1}{2}\begin{bmatrix}x\\y\end{bmatrix}^T\begin{bmatrix}H_1 & H_2\\H_2^T & H_3\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix} + \begin{bmatrix}e_1^T&e_2^T\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}\\ & \text{subject to } \, \begin{bmatrix}F_1 & F_2\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}\leq h \end{aligned} \end{aligned}\]The inner problem constraining the follower \(y\), is limited to convex quadratic programming problems. The outer problem is allowed to be essentially anything that YALMIP can handle.
The bilevel solver that is available in YALMIP replaces the optimality condition on \(y^{\star}\) with the KKT conditions.
\[\begin{aligned} H_3y + H_2^Tx + e_2 - F_2^T\lambda &=0\\ \lambda &\geq 0\\ h-F_1x-F_2y & \geq 0\\ \lambda^T(h-F_1x-F_2y) & = 0 \end{aligned}\]This is precisely what is done in the manually derived bilevel solution methods in bilevel example, but the possible benefit of using YALMIPs native support as we will do here is that this solver branches directly on the complementarity conditions, and thus avoids to introduce any numerically dangerous big-M constants.
Bilevel linear and quadratic programming
Let us start with a simple bilevel linear programming problem. Start by defining the outer (leader) variables \(x\) and inner (follower) variables \(y\).
sdpvar x1 x2
sdpvar y1 y2 y3
Using standard notation, we define the outer constraints and objective
OO = -8*x1 - 4*x2 + 4*y1 - 40*y2 + 4*y3;
CO = [[y1 y2 y3]>=0,[x1 x2]>=0];
followed by the inner constraints and objective
OI = x1 + 2*x2 + y1 + y2 + 2*y3;
CI = [-y1 + y2 + y3 <= 1,
2*x1 - y1 + 2*y2 - 0.5*y3 <= 1,
2*x2 + 2*y1 - y2 - 0.5*y3 <= 1,
[y1 y2 y3] >= 0,
[x1 x2] >= 0];
Having defined the optimization structures, we call the bilevel solver with the constraints and objectives, and the variables \(y\) in order to tell the solver which variables are the inner variables.
solvebilevel(CO,OO,CI,OI,[y1 y2 y3]);
* Starting YALMIP bilevel solver.
* Outer solver : GUROBI-GUROBI
* Inner solver : GUROBI-GUROBI
* Max iterations : 300
Node Upper Gap(%) Lower Open
1 : Inf Inf -5.000E+01 2 Solved to optimality
2 : Inf Inf -5.000E+01 3 Solved to optimality
3 : -6.000E+00 78.57 -5.000E+01 2 Solved to optimality
4 : -6.000E+00 78.57 -5.000E+01 1 Infeasible
5 : -6.000E+00 62.50 -2.600E+01 2 Solved to optimality
6 : -2.600E+01 0.00 -2.600E+01 0 Solved to optimality
As you hopefully have noticed, completely standard YALMIP code is used to setup and manipulate the model. Hence, to obtain the final solution, we use value.
value([y1 y2 y3])
ans =
0 0.6000 0.4000
Adding additional complicating constraints is allowed, as long as YALMIP can identify a solver which is capable of solving the outer problem appended with the KKT conditions, excluding the nonconvex complementary slackness constraint. Hence, we can add integrality constraints easily to our model
solvebilevel([CO, integer(y2)],OO,CI,OI,[y1 y2 y3]);
In the same sense, convex quadratic problems are dealt with straightforwardly
solvebilevel(CO,OO+OO^2,CI,OI^2,[y1 y2 y3]);
Bilevel programming with general outer problem
A strong feature of the built-in solver is that it builds upon the infrastructure in YALMIP, and easily hooks up to almost any kind of outer problem. Hence, we can take the problem above, and append a semidefinite constraint to the outer problem. The only difference is that during the branching of the complementary conditions, semidefinite programs have to be solved in each node.
CO = [[y1 y2 y3] >= 0,[x1 x2] >= 0];
CO = [CO, [1 x1+x2;x1+x2 1/2] >= 0];
solvebilevel(CO,OO,CI,OI,[y1 y2 y3]);
* Starting YALMIP bilevel solver.
* Outer solver : SeDuMi-1.1
* Inner solver : GLPK-GLPKMEX
* Max iterations : 300
Node Upper Gap(%) Lower Open
1 : 3.254E-015 100.00 -5.000E+001 2 Solved to optimality
2 : -1.723E-010 100.00 -5.000E+001 3 Solved to optimality
3 : -4.828E+000 82.39 -5.000E+001 2 Solved to optimality
4 : -4.828E+000 82.39 -5.000E+001 1 Infeasible in solver
5 : -1.940E+001 13.07 -2.523E+001 2 Solved to optimality
6 : -1.940E+001 13.07 -2.523E+001 3 Solved to optimality
7 : -1.940E+001 13.07 -2.523E+001 4 Solved to optimality
8 : -1.940E+001 13.07 -2.523E+001 3 Infeasible in solver
9 : -1.940E+001 7.18 -2.240E+001 4 Solved to optimality
10 : -1.940E+001 7.18 -2.240E+001 3 Infeasible in solver
11 : -1.940E+001 0.00 -1.940E+001 2 Infeasible in solver
The bilevel solver is restricted to convex quadratic inner problems, but convexity is not a requirement on the outer problems. The bilevel solver solves the outer problem repeatedly in a branch-and-bound procedure, with additional equality constraints derived from complementary slackness appended. Hence, we have to make a choice between global solution of the outer problem, or a simple local solution. If we go for a local solution (using, e.g, fmincon), we have no guarantees (except that the inner optimality constraint is satisfied). If we use a global outer solver (such as bmibnb), a globally optimal bilevel solution follows.
As an illustration, we solve the original problem, but append a nonconvex quadratic term to the outer problem. To ensure a globally optimal solution, we use the global solver bmibnb as the outer solver.
OO = -8*x1 - 4*x2 + 4*y1 - 40*y2 + 4*y3;
CO = [[y1 y2 y3] >= 0,[x1 x2] >= 0];
OI = x1 + 2*x2 + y1 + y2 + 2*y3;
CI = [-y1 + y2 + y3 <= 1,
2*x1 - y1 + 2*y2 - 0.5*y3 <= 1,
2*x2 + 2*y1 - y2 - 0.5*y3 <= 1,
[y1 y2 y3] >= 0,
[x1 x2] >= 0];
ops = sdpsettings('bilevel.outersolver','bmibnb');
solvebilevel(CO,OO-OO^2,CI,OI,[y1 y2 y3],ops)
* Starting YALMIP bilevel solver.
* Outer solver : BMIBNB
* Inner solver : GLPK-GLPKMEX
* Max iterations : 300
Node Upper Gap(%) Lower Open
1 : -1.110E-015 100.00 -2.550E+003 2 Solved to optimality
2 : -1.110E-015 100.00 -2.550E+003 3 Solved to optimality
3 : -1.110E-015 100.00 -2.550E+003 4 Solved to optimality
4 : -1.110E-015 100.00 -2.550E+003 5 Solved to optimality
5 : -7.020E+002 56.83 -2.550E+003 4 Solved to optimality
6 : -7.020E+002 56.83 -2.550E+003 3 Infeasible in solver
7 : -7.020E+002 30.97 -1.332E+003 2 Solved to optimality
8 : -7.020E+002 30.97 -1.332E+003 3 Solved to optimality
9 : -7.020E+002 30.97 -1.332E+003 2 Solved to optimality
10 : -7.020E+002 30.97 -1.332E+003 1 Infeasible in solver
11 : -7.020E+002 0.00 -7.020E+002 0 Solved to optimality
