Added 18/06/2025
Archetype
linear_bilevel
Datasets
Description
Anandalingam, G. and White, D.J. (1990) (see page 1172).
A solution method for the linear static Stackelberg problem using penalty functions.
https://doi.org/10.1109/9.58565Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 6,
"h": 0
}Solution
{
"optimality": "global",
"x": [16],
"y": [11],
"F": -49,
"G": [16],
"H": [],
"f": 17,
"g": [28,12,0,0,12,11],
"h": []
}Description
Bard, Jonathan F. (1984) (see page 18).
Optimality conditions for the bilevel programming problem.
https://onlinelibrary.wiley.com/doi/abs/10.1002/nav.3800310104Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 5,
"h": 0
}Solution
{
"optimality": "global",
"x": [0.8888888888888888],
"y": [2.2222222222222223],
"F": 3.111111111111111,
"G": [0.8888888888888888],
"H": [],
"f": -6.666666666666667,
"g": [0,0,6,7.555555555555555,2.2222222222222223],
"h": []
}Description
Bard, Jonathan F. (1984) (see page 18).
Optimality conditions for the bilevel programming problem.
https://onlinelibrary.wiley.com/doi/abs/10.1002/nav.3800310104Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 5,
"h": 0
}Solution
{
"optimality": "best known",
"x": [7.2],
"y": [1.6],
"F": -37.6,
"G": [7.2],
"H": [],
"f": 1.6,
"g": [6,2.2,0,0,1.6],
"h": []
}Description
Bard, Jonathan F. (1991) (see page 374).
Some properties of the bilevel programming problem.
https://doi.org/10.1007/BF00941574Dimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 5,
"h": 0
}Solution
{
"optimality": "global",
"x": [0],
"y": [0,1],
"F": -1,
"G": [0],
"H": [],
"f": -1,
"g": [1,0,0,0,1],
"h": []
}Description
Bard, Jonathan F. and Falk, James E. (1982) (see page 90).
An explicit solution to the multi-level programming problem.
https://doi.org/10.1016/0305-0548(82)90007-7Dimension
{
"x": 2,
"y": 2,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 5,
"h": 0
}Solution
{
"optimality": "global",
"x": [2,0],
"y": [1.5,0],
"F": -3.25,
"G": [2,0],
"H": [],
"f": -4,
"g": [0,0,0,1.5,0],
"h": []
}Description
Ben-Ayed, Omar and Blair, E. (1990) (see page 557).
Computational Difficulties of Bilevel Linear Programming.
https:/doi.org/10.1287/opre.38.3.556Dimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 4,
"h": 0
}Solution
{
"optimality": "global",
"x": [1],
"y": [0,1],
"F": -2.5,
"G": [1,0],
"H": [],
"f": -5,
"g": [3,0,0,1],
"h": []
}Description
Ben-Ayed, Omar and Blair, E. (1990) (see page 558).
Computational Difficulties of Bilevel Linear Programming.
https:/doi.org/10.1287/opre.38.3.556Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 4,
"h": 0
}Solution
{
"optimality": "global",
"x": [1],
"y": [5],
"F": -6,
"G": [1],
"H": [],
"f": 5,
"g": [0,0,5,5],
"h": []
}Description
Bialas, Wayne F and Karwan, Mark H (1984) (see page 1009).
Two-level linear programming.
https://www.jstor.org/stable/2631591Dimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 7,
"h": 0
}Solution
{
"optimality": "global",
"x": [1.5],
"y": [1,0.5],
"F": -2,
"G": [1.5],
"H": [],
"f": -0.5,
"g": [0,1,1,0,0,1,0.5],
"h": []
}Description
Bialas, Wayne F and Karwan, Mark H (1984) (see page 1016).
Two-level linear programming.
https://www.jstor.org/stable/2631591Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 6,
"h": 0
}Solution
{
"optimality": "global",
"x": [16],
"y": [11],
"F": -11,
"G": [16],
"H": [],
"f": 11,
"g": [28,0,12,12,0,11],
"h": []
}Description
Candler, Wilfred and Townsley, Robert (1982) (see page 91).
A linear two-level programming problem.
https://doi.org/10.1016/0305-0548(82)90006-5Dimension
{
"x": 2,
"y": 3,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 6,
"h": 0
}Solution
{
"optimality": "global",
"x": [0,0.9],
"y": [0,0.6,0.4],
"F": -29.2,
"G": [0,0.9],
"H": [],
"f": 3.2,
"g": [0,0,0,0,0.6,0.4],
"h": []
}Description
Clark, PA and Westerberg, AW (1988) (see page 414).
A note on the optimality conditions for the bilevel programming problem.
https://doi.org/10.1002/1520-6750(198810)35:5%3C413::AID-NAV3220350505%3E3.0.CO;2-6Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 0,
"H": 0,
"f": 1,
"g": 3,
"h": 0
}Solution
{
"optimality": "global",
"x": [19],
"y": [14],
"F": -37,
"G": [],
"H": [],
"f": 14,
"g": [24,0,0],
"h": []
}Description
Clark, Peter A and Westerberg, Arthur W (1990) (see page 89).
Bilevel programming for steady-state chemical process design-I. Fundamentals and algorithms.
https://doi.org/10.1016/0098-1354(90)87007-CDimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 5,
"h": 0
}Solution
{
"optimality": "best known",
"x": [5],
"y": [4,2],
"F": -13,
"G": [5,3],
"H": [],
"f": -4,
"g": [4,0,14,0,0],
"h": []
}Description
Glackin, J and Ecker, JG and Kupferschmid, M (2009) (see page 206).
Solving bilevel linear programs using multiple objective linear programming.
https://doi.org/10.1007/s10957-008-9467-2Dimension
{
"x": 2,
"y": 1,
"F": 1,
"G": 3,
"H": 0,
"f": 1,
"g": 3,
"h": 0
}Solution
{
"optimality": "global",
"x": [1,2],
"y": [0],
"G": [0,1,2],
"H": [],
"g": [1,0,0],
"h": []
}Description
Haurie, Alain and Savard, G and White, Douglas J (1990) (see page 554).
A note on: an efficient point algorithm for a linear two-stage optimization problem.
https://doi.org/10.1287/opre.38.3.553Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 0,
"H": 0,
"f": 1,
"g": 4,
"h": 0
}Solution
{
"optimality": "global",
"x": [12],
"y": [3],
"F": 27,
"G": [],
"H": [],
"f": -3,
"g": [36,0,0,7],
"h": []
}Description
Hu, Tiesong and Guo, Xuning and Fu, Xiang and Lv, Yibing (2010) (see page 241).
A neural network approach for solving linear bilevel programming problem.
https://doi.org/10.1016/j.knosys.2010.01.001Dimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 5,
"h": 0
}Solution
{
"optimality": "global",
"x": [1.8888888888888888],
"y": [0.8888888888888888,0],
"F": -8.444444444444445,
"G": [1.8888888888888888],
"H": [],
"f": -4.555555555555555,
"g": [0,0.11111111111111116,0,0.8888888888888888,0],
"h": []
}Dimension
{
"x": 10,
"y": 10,
"F": 1,
"G": 0,
"H": 0,
"f": 1,
"g": 11,
"h": 0
}Solution
{
"optimality": "global",
"x": [1,1,1,1,1,1,1,1,1,1],
"y": [1,1,1,1,1,1,1,1,1,1],
"F": 10,
"G": [],
"H": [],
"f": 10,
"h": []
}Dimension
{
"x": 100,
"y": 100,
"F": 1,
"G": 0,
"H": 0,
"f": 1,
"g": 101,
"h": 0
}Solution
{
"optimality": "global",
"F": 100,
"G": [],
"H": [],
"f": 100,
"h": []
}Dimension
{
"x": 3,
"y": 3,
"F": 1,
"G": 0,
"H": 0,
"f": 1,
"g": 4,
"h": 0
}Solution
{
"optimality": "global",
"x": [1,1,1],
"y": [1,1,1],
"F": 3,
"G": [],
"H": [],
"f": 3,
"g": [0,0,0,0],
"h": []
}Dimension
{
"x": 50,
"y": 50,
"F": 1,
"G": 0,
"H": 0,
"f": 1,
"g": 51,
"h": 0
}Solution
{
"optimality": "global",
"F": 50,
"G": [],
"H": [],
"f": 50,
"h": []
}Description
Lan, Kuen-Ming and Wen, Ue-Pyng and Shih, Hsu-Shih and Lee, E Stanley (2007) (see page 882).
A hybrid neural network approach to bilevel programming problems.
https://doi.org/10.1016/j.aml.2006.07.013Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 7,
"h": 0
}Solution
{
"optimality": "best known",
"x": [17.45],
"y": [10.908],
"F": -85.088,
"G": [17.45],
"H": [],
"f": 50.174,
"g": [8.366,0.008000000000002672,0.018000000000000682,32.194,80.25999999999999,53.081999999999994,10.908],
"h": []
}Description
Liu, Yi-Hsin and Hart, Stephen M (1994) (see page 166).
Characterizing an optimal solution to the linear bilevel programming problem.
https://doi.org/10.1016/0377-2217(94)90155-4Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 4,
"h": 0
}Solution
{
"optimality": "global",
"x": [4],
"y": [4],
"F": -16,
"G": [4],
"H": [],
"f": 4,
"g": [3,0,0,4],
"h": []
}Description
Mersha, Ayalew Getachew and Dempe, Stephan (2006) (see page 250).
Linear bilevel programming with upper level constraints depending on the lower level solution.
https://doi.org/10.1016/j.amc.2005.11.134Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 5,
"h": 0
}Solution
{
"optimality": "best known",
"x": [9],
"y": [6],
"F": -39,
"G": [9],
"H": [],
"f": 6,
"g": [0,12,50,0,6],
"h": []
}Description
Mersha, Ayalew Getachew and Dempe, Stephan (2006) (see page 251).
Linear bilevel programming with upper level constraints depending on the lower level solution.
https://doi.org/10.1016/j.amc.2005.11.134Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [8],
"y": [6],
"F": -20,
"G": [10,0],
"H": [],
"f": -6,
"g": [15,0],
"h": []
}Description
Tuy, Hoang and Migdalas, Athanasios and Värbrand, Peter (1993) (see page 17).
A global optimization approach for the linear two-level program.
https://doi.org/10.1007/BF01100237Dimension
{
"x": 2,
"y": 2,
"F": 1,
"G": 3,
"H": 0,
"f": 1,
"g": 4,
"h": 0
}Solution
{
"optimality": "global",
"x": [2,0],
"y": [1.5,0],
"F": -3.25,
"G": [0,2,0],
"H": [],
"f": -6,
"g": [0,0,1.5,0],
"h": []
}Description
Tuy, Hoang and Migdalas, Athanasios and Värbrand, Peter (1994) (see page 257).
A quasiconcave minimization method for solving linear two-level programs.
https://doi.org/10.1007/BF01098360Dimension
{
"x": 2,
"y": 2,
"F": 1,
"G": 3,
"H": 0,
"f": 1,
"g": 3,
"h": 0
}Solution
{
"optimality": "global",
"x": [0,3],
"y": [0,0],
"F": 6,
"G": [1,0,3],
"H": [],
"f": 0,
"g": [0,0,0],
"h": []
}Description
Tuy, Hoang and Migdalas, Athanasios and Värbrand, Peter (1993) (see page 551).
A global optimization approach for the linear two-level program.
https://doi.org/10.1007/BF01100237Dimension
{
"x": 10,
"y": 6,
"F": 1,
"G": 22,
"H": 0,
"f": 1,
"g": 19,
"h": 0
}Solution
{
"optimality": "infeasible",
"x": [0,8.170692,10,0,7.27894,3.042311,0,10,0.001982,9.989153],
"y": [3.10128,10,10,10,0,9.846133],
"F": -467.461261,
"H": [],
"f": -11.619362,
"h": []
}Description
Visweswaran, V and Floudas, CA and Ierapetritou, MG and Pistikopoulos, EN (1996) (see page 159).
A decomposition-based global optimization approach for solving bilevel linear and quadratic programs.
https://doi.org/10.1007/978-1-4613-3437-8_10Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 5,
"h": 0
}Solution
{
"optimality": "best known",
"x": [0.8888888888888888],
"y": [2.2222222222222223],
"F": 3.111111111111,
"G": [0.8888888888888888],
"H": [],
"f": -6.666666666666,
"g": [0,0,6,5.555555555555555,2.2222222222222223],
"h": []
}Description
Wang, Yuping and Jiao, Yong-Chang and Li, Hong (2005) (see page 228).
An evolutionary algorithm for solving nonlinear bilevel programming based on a new constraint-handling scheme.
https://doi.org/10.1109/TSMCC.2004.841908Dimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "best known",
"x": [0],
"y": [1,0],
"F": -1000,
"G": [0,1],
"H": [],
"f": -1,
"g": [0,0],
"h": []
}Description
Thürauf, Johannes and Kleinert, Thomas and Ljubić, Ivana and Ralphs, Ted and Schmidt, Martin (2024).
BOBILib: Bilevel Optimization (Benchmark) Instance Library.
https://optimization-online.org/?p=27063Dimension
{
"x": 19,
"y": 8,
"F": 1,
"G": 39,
"H": 0,
"f": 1,
"g": 44,
"h": 0
}Solution
{
"optimality": "unknown"
}Description
Thürauf, Johannes and Kleinert, Thomas and Ljubić, Ivana and Ralphs, Ted and Schmidt, Martin (2024).
BOBILib: Bilevel Optimization (Benchmark) Instance Library.
https://optimization-online.org/?p=27063Dimension
{
"x": 40,
"y": 40,
"F": 1,
"G": 100,
"H": 0,
"f": 1,
"g": 110,
"h": 0
}Solution
{
"optimality": "unknown"
}$title linear_bilevel
$onText
Upper-level
============
minimise c_upper ' x + d_upper ' y
subject to A_upper * x + B_upper * y <= rhs_upper
x >= x_lower_bound
x <= x_upper_bound
y solves lower-level
Lower-level
============
minimise c_lower ' x + d_lower ' y
subject to A_lower * x + B_lower * y <= rhs_lower
y >= y_lower_bound
y <= y_upper_bound
Table of Dimensions
===================
| Array | Shape |
| --------- | -------------- |
| A_upper | (n_upper, n_x) |
| B_upper | (n_upper, n_y) |
| c_upper | (n_x,) |
| d_upper | (n_y,) |
| rhs_upper | (n_upper,) |
| A_lower | (n_lower, n_x) |
| B_lower | (n_lower, n_y) |
| c_lower | (n_x,) |
| d_lower | (n_y,) |
| rhs_lower | (n_lower,) |
$offText
* Sets
Sets
i "Index for x variables"
j "Index for y variables"
m "Index for upper-level inequality constraints"
l "Index for Lower-level inequality constraints";
* Parameters
Parameters
A_upper(m,i) "Upper-level inequality matrix A"
B_upper(m,j) "Upper-level inequality matrix B"
c_upper(i) "Upper-level objective coefficients for x"
d_upper(j) "Upper-level objective coefficients for y"
rhs_upper(m) "Upper-level inequality right hand side"
A_lower(l,i) "Lower-level inequality matrix A"
B_lower(l,j) "Lower-level inequality matrix B"
c_lower(i) "Lower-level objective coefficients for x"
d_lower(j) "Lower-level objective coefficients for y"
rhs_lower(l) "Lower-level inequality right hand side"
x_lower_bound(i) "x lower bound"
x_upper_bound(i) "x upper bound"
y_lower_bound(j) "y lower bound"
y_upper_bound(j) "y upper bound";
* Variables
Variables
x(i) "Upper-level decision variables"
y(j) "Lower-level decision variables"
obj "Objective function value";
* Data
$gdxin "/bolib3/data/linear_bilevel_gdx/tuy_et_al_2007_ex3.gdx"
$load i=i, j=j, m=m, l=l,
$load A_upper=A_upper, B_upper=B_upper, c_upper=c_upper, d_upper=d_upper, rhs_upper=rhs_upper, x_lower_bound, x_upper_bound
$load A_lower=A_lower, B_lower=B_lower, c_lower=c_lower, d_lower=d_lower, rhs_lower=rhs_lower, y_lower_bound, y_upper_bound
Display c_upper;
Display c_upper;
Variables x, y, obj_val_upper, obj_val_lower;
* Equations
Equations
obj_eq_upper "Upper-level objective"
obj_eq_lower "Lower-level objective"
upper_ineq(m) "Upper-level inequality constraints"
lower_ineq(l) "Lower-level inequality constraint";
* Upper-level objective
obj_eq_upper .. obj_val_upper =e= sum(i, c_upper(i)*x(i)) + sum(j, d_upper(j)*y(j));
* Lower-level objective
obj_eq_lower .. obj_val_lower =e= sum(i, c_lower(i)*x(i)) + sum(j, d_lower(j)*y(j));
* Upper-level inequality constraints
upper_ineq(m) .. sum(i, A_upper(m,i)*x(i)) + sum(j, B_upper(m,j)*y(j)) =g= rhs_upper(m);
* Lower-level inequality constraint
lower_ineq(l) .. sum(i, A_lower(l,i)*x(i)) + sum(j, B_lower(l,j)*y(j)) =g= rhs_lower(l);
* Bounds
x.lo(i) = x_lower_bound(i);
x.up(i) = x_upper_bound(i);
y.lo(j) = y_lower_bound(j);
y.up(j) = y_upper_bound(j);
* Solve
model linear_bilevel / all /;
$echo bilevel x min obj_val_lower y obj_eq_lower, lower_ineq > "%emp.info%"
solve linear_bilevel us emp min obj_val_upper;
\subsection{linear\_bilevel}
\label{subsec:linear_bilevel}
% Linear Bilevel Program
For decision variables $\x\in\R^{n_x}$ and $\y\in\R^{n_y}$ the \emph{(fully) Linear Bilevel Program}~\eqref{eq:LBP} is defined as:
\begin{flalign*}
\label{eq:LBP}
\tag{LBP}
\minimise_{\x\in\R^{n_x}, \y^*\in\R^{n_y}} \quad
& {c^\top_\upper} \x + {d^\top_\upper} \y^*\\
\subjectto \quad
& {A_\upper} \x + {B_\upper} \y^* \geq {r_\upper},\\
& \xlb\leq \x \leq\xub,\\
& y^* \in \argmin_{y\in\R^{n_y}}
\left\{
\begin{array}{l}
{c^\top_\low} \x + {d^\top_\low} \y\\
\st[1]\quad {A_\low} \x + {B_\low} \y \geq {r_\low},\\
\st[2]\quad \ylb\leq \y \leq\yub.
\end{array}
\right.
\end{flalign*}
where the program is parametrised by the coefficients listed in Table~\ref{tab:parameters_LBP}.
% //==================================================\\
% || Table: Parameters Linear Bilevel ||
% \\==================================================//
\begin{table}[H]
\centering
\begin{tabular}{c c l}
\toprule
Parameter & Varaible name & Dimension \\
\midrule
$A_\upper$ & \texttt{"A\_upper"} & $\in\R^{n_G \times n_x}$ \\
$B_\upper$ & \texttt{"B\_upper"} & $\in\R^{n_G \times n_y}$ \\
$c_\upper$ & \texttt{"c\_upper"} & $\in\R^{n_x}$ \\
$d_\upper$ & \texttt{"d\_upper"} & $\in\R^{n_y}$ \\
$r_\upper$ & \texttt{"rhs\_upper"} & $\in\R^{n_G}$ \\
\midrule
$A_\low$ & \texttt{"A\_lower"} & $\in\R^{n_g \times n_x}$ \\
$B_\low$ & \texttt{"B\_lower"} & $\in\R^{n_g \times n_y}$ \\
$c_\low$ & \texttt{"c\_lower"} & $\in\R^{n_x}$ \\
$d_\low$ & \texttt{"d\_lower"} & $\in\R^{n_y}$ \\
$r_\low$ & \texttt{"rhs\_lower"} & $\in\R^{n_g}$ \\
\midrule
$\xlb$ & \texttt{"x\_lower\_bound"} & $\in\R^{n_x}$\\
$\xub$ & \texttt{"x\_upper\_bound"} & $\in\R^{n_x}$\\
$\ylb$ & \texttt{"y\_lower\_bound"} & $\in\R^{n_y}$\\
$\yub$ & \texttt{"y\_upper\_bound"} & $\in\R^{n_y}$\\
\bottomrule
\end{tabular}
\caption{Parameters of the \eqref{eq:LBP}.}
\label{tab:parameters_LBP}
\end{table}
% //==================================================\\
% || Data Files ||
% \\==================================================//
We provide data files containing the parameters for over 60 instances of~\eqref{eq:LBP}.
At the time of writing, the library supports two different formats:
\begin{itemize}
\itemsep0em
% JSON
\item \textbf{JSON}.
This is a key-value format where the key is one of the variable names given in Table~\ref{tab:parameters_LBP} and the value is a vector/matrix of coefficients.
See \href{https://github.com/bolib3/bolib3/blob/main/automation/templates/archetype_templates/linear_template.json}{linear\_template.json}.
% AUX + MPS
\item \textbf{AUX + MPS}.
This is a pairing of a Mathematical Programming System (MPS) file for storing and recording linear constraints and an Auxiliary (AUX) file for signalling which variables and constraints are lower-level.
It was developed for the COIN-OR Mibs project~\cite[page 552, Design of Mibs]{Tahernejad2020}.
For an in-depth documentation of this file format, read \href{https://coin-or.github.io/MibS/input.html}{coin-or.github.io/MibS/input.html}.
\end{itemize}classdef linear_bilevel
% Linear Bilevel Program (LBP)
% ============================
% The LBP, sometimes called fully-linear, has both linear objective functions and linear constraints.
% Each instance is defined by vectors of coefficients stored in bolib3/data/linear_bilevel.
%
%
% Upper-level
% ===========
% minimise c_upper ' x + d_upper ' y
% subject to A_upper * x + B_upper * y >= rhs_upper
% x >= x_lower_bound
% x <= x_upper_bound
% y solves lower-level
%
%
% Lower-level
% ===========
% minimise c_lower ' x + d_lower ' y
% subject to A_lower * x + B_lower * y >= rhs_lower
% y >= y_lower_bound
% y <= y_upper_bound
%
%
% Table of Dimensions
% ===================
% | Array | Shape |
% | --------- | ---------- |
% | A_upper | (n_G, n_x) |
% | B_upper | (n_G, n_y) |
% | c_upper | (n_x,) |
% | d_upper | (n_y,) |
% | rhs_upper | (n_G,) |
% | A_lower | (n_g, n_x) |
% | B_lower | (n_g, n_y) |
% | c_lower | (n_x,) |
% | d_lower | (n_y,) |
% | rhs_lower | (n_g,) |
% were n_x and n_y are the number of upper-level and lower-level decision variables respectively.
% and n_G and n_g are the number of upper-level and lower-level inequality constraints respectively.
%
%
% Example MATLAB use case
% =======================
% addpath('bolib3/matlab');
% data = linear_bilevel.read_data('bolib3/data/linear_bilevel/anandalingham_white_1990.json');
% x = [16];
% y = [11];
% linear_bilevel.g(x, y, data)
%
%
% Accepted data formats
% =====================
% * JSON: A key-value format where the key is a variable name and the value is a vector of coefficients.
% See the function _read_json() for more details.
% * AUX + MPS: The auxiliary (AUX) file lists names of the lower-level variables and constraints.
% The Mathematical Programming System (MPS) file defines all constraints in a name-based system.
% See the function _read_aux_mps() for more details.
%
properties(Constant)
name = 'linear_bilevel';
category = 'archetype';
subcategory = '';
author = 'Samuel Ward';
datasets = sort({
'anandalingham_white_1990.json'; 'bard_1984a.json'; 'bard_1984b.json'; 'bard_1991_ex2.json';
'bard_falk_1982_ex2.json'; 'ben_ayed_blair_1990a.json'; 'ben_ayed_blair_1990b.json';
'bialas_karwan_1984a.json'; 'bialas_karwan_1984b.json'; 'candler_townsley_1982.json';
'clark_westerberg_1988.json'; 'clark_westerberg_1990b.json'; 'glackin_et_al_2009.json';
'haurie_savard_white_1990.json'; 'hu_huang_zhang_2009.json'; 'hypercube_10.json'; 'hypercube_100.json';
'hypercube_3.json'; 'hypercube_50.json'; 'lan_wen_shih_lee_2007.json'; 'liu_hart_1994.json';
'mersha_dempe_2006_ex1.json'; 'mersha_dempe_2006_ex2.json';'tuy_et_al_1993.json'; 'tuy_et_al_1994.json';
'tuy_et_al_2007_ex3.json'; 'visweswaran_et_al_1996.json'; 'wang_jiao_li_2005.json'
});
paths = fullfile('bolib3', 'data', 'linear_bilevel', linear_bilevel.datasets);
end
methods(Static)
function obj = F(x, y, data)
% Upper-level objective function
% minimise F(x,y) := c'x + d'y
obj = dot(data.c_upper, x) + dot(data.d_upper, y);
end
function ineq = G(x, y, data)
% Upper-level inequality constraints
% Ax + By >= rhs
% x >= x_lower_bound
% x <= x_upper_bound
ineq = [];
if ~isempty(data.rhs_upper)
ineq = [ineq; data.A_upper * x + data.B_upper * y - data.rhs_upper];
end
if ~isempty(data.x_lower_bound)
ineq = [ineq; x - data.x_lower_bound];
end
if ~isempty(data.x_upper_bound)
ineq = [ineq; data.x_upper_bound - x];
end
end
function val = H(~, ~, ~)
% Upper-level equality constraints
val = [];
end
function obj = f(x, y, data)
% Lower-level objective function
% minimise f(x,y) := c'x + d'y
obj = dot(data.c_lower, x) + dot(data.d_lower, y);
end
function ineq = g(x, y, data)
% Lower-level inequality constraints
ineq = [];
if ~isempty(data.rhs_lower)
ineq = [ineq; data.A_lower * x + data.B_lower * y - data.rhs_lower];
end
if ~isempty(data.y_lower_bound)
ineq = [ineq; y - data.y_lower_bound];
end
if ~isempty(data.y_upper_bound)
ineq = [ineq; data.y_upper_bound - y];
end
end
function val = h(~, ~, ~)
% Lower-level equality constraints
val = [];
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=- %
%% Derivatives %%
% -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=- %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
function dx = dF_dx(~, ~, data)
% ∂F/∂x
% Dimention (n_x × 1)
dx = data.c_upper;
end
function dy = dF_dy(~, ~, data)
% ∂F/∂y
% Dimention (n_y × 1)
dy = data.d_upper;
end
function Hxx = d2F_dxx(~, ~, data)
% ∂²F/∂x²
% Dimention (n_x × n_x)
Hxx = zeros(data.n_x, data.n_x);
end
function Hyy = d2F_dyy(~, ~, data)
% ∂²F/∂y²
% Dimention (n_y × n_y)
Hyy = zeros(data.n_y, data.n_y);
end
function Hxy = d2F_dxy(~, ~, data)
% ∂²F/∂x∂y
% Dimention (n_x × n_y)
Hxy = zeros(data.n_x, data.n_y);
end
function Jx = dG_dx(~, ~, data)
% Jacobian of G wrt x
% Dimention (|G| × n_x)
Jx = zeros(0, data.n_x);
if ~isempty(data.rhs_upper)
Jx = [Jx; data.A_upper];
end
if ~isempty(data.x_lower_bound)
Jx = [Jx; eye(data.n_x)];
end
if ~isempty(data.x_upper_bound)
Jx = [Jx; -eye(data.n_x)];
end
end
function Jy = dG_dy(~, ~, data)
% Jacobian of G wrt y:
% Dimention (|G| × n_y)
Jy = zeros(0, data.n_y);
if ~isempty(data.rhs_upper)
Jy = [Jy; data.B_upper];
end
if ~isempty(data.x_lower_bound)
Jy = [Jy; zeros(data.n_x, data.n_y)];
end
if ~isempty(data.x_upper_bound)
Jy = [Jy; zeros(data.n_x, data.n_y)];
end
end
function Jx = dH_dx(~, ~, data)
% Jacobian of H wrt x
% Dimention (0 × n_x)
Jx = zeros(0, data.n_x);
end
function Jy = dH_dy(~, ~, data)
% Jacobian of H wrt y
% Dimention (0 × n_y)
Jy = zeros(0, data.n_y);
end
function dx = df_dx(~, ~, data)
% ∂f/∂x
% Dimention (n_x × 1)
dx = data.c_lower;
end
function dy = df_dy(~, ~, data)
% ∂f/∂y
% Dimention (n_y × 1)
dy = data.d_lower;
end
function Hxx = d2f_dxx(~, ~, data)
% ∂²f/∂x²
% Dimention (n_x × n_x)
Hxx = zeros(data.n_x, data.n_x);
end
function Hyy = d2f_dyy(~, ~, data)
% ∂²f/∂y²
% Dimention (n_y × n_y)
Hyy = zeros(data.n_y, data.n_y);
end
function Hxy = d2f_dxy(~, ~, data)
% ∂²f/∂x∂y
% Dimention (n_x × n_y)
Hxy = zeros(data.n_x, data.n_y);
end
function Jx = dg_dx(~, ~, data)
% Jacobian of g wrt x
% Dimention (|g| × n_x)
Jx = zeros(0, data.n_x);
if ~isempty(data.rhs_lower)
Jx = [Jx; data.A_lower];
end
if ~isempty(data.y_lower_bound)
Jx = [Jx; zeros(data.n_y, data.n_x)];
end
if ~isempty(data.y_upper_bound)
Jx = [Jx; zeros(data.n_y, data.n_x)];
end
end
function Jy = dg_dy(~, ~, data)
% Jacobian of g wrt y
% Dimention (|g| × n_y)
Jy = zeros(0, data.n_y);
if ~isempty(data.rhs_lower)
Jy = [Jy; data.B_lower];
end
if ~isempty(data.y_lower_bound)
Jy = [Jy; eye(data.n_y)];
end
if ~isempty(data.y_upper_bound)
Jy = [Jy; -eye(data.n_y)];
end
end
function Jx = dh_dx(~, ~, data)
% Jacobian of h wrt x
% Dimention (0 × n_x)
Jx = zeros(0, data.n_x);
end
function Jy = dh_dy(~, ~, data)
% Jacobian of h wrt y
% Dimention (0 × n_y)
Jy = zeros(0, data.n_y);
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=- %
%% Read Data %%
% -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=- %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
function data = read_data(path)
% READ_DATA
% data = linear_bilevel.read_data('.../example.json') % JSON
% data = linear_bilevel.read_data('.../example.aux') % AUX ... will find example.mps or example.mps.gz
% JSON
if endsWith(path, ".json", IgnoreCase=true)
data = linear_bilevel.read_json(path);
% AUX
% Find the matching .mps or .mps.gz in the same folder
elseif endsWith(path, ".aux", IgnoreCase=true)
mps_path = extractBefore(string(path), ".aux") + ".mps";
mps_gz_path = extractBefore(string(path), ".aux") + ".mps.gz";
if exist(mps_path, 'file')
data = linear_bilevel.read_aux_mps(path, which(mps_path));
elseif exist(mps_gz_path, 'file')
data = linear_bilevel.read_aux_mps(path, which(mps_gz_path));
else
error( ...
"Failed to find" + ...
sprintf("\n (.mps): ") + mps_path + ...
sprintf("\n or (.mps.gz): ") + mps_gz_path + ...
sprintf("\n Working directory: ") + string(pwd) ...
);
end
% OTHERWISE
else
error("Path must be .json or .aux: got '%s'.", path);
end
% Final dimension check
linear_bilevel.validate_data_dimensions(data);
end
function validate_data_dimensions(data)
% VALIDATE_DATA_DIMENSIONS
% See the table of dimensions in the docstring at the top of this file.
n_x = data.n_x;
n_y = data.n_y;
n_G = numel(data.rhs_upper);
n_g = numel(data.rhs_lower);
mustBeSize = @(X, expected_size, name, allow_empty) assert( ...
isequal(size(X), expected_size) || (isempty(X) && allow_empty), ...
"Dimension of %s should be %s but got %s. See: '%s'.", ...
name, mat2str(expected_size), mat2str(size(X)), string(data.path) ...
);
% Upper-level arrays
mustBeSize(data.A_upper, [n_G, n_x], 'A_upper', 0);
mustBeSize(data.B_upper, [n_G, n_y], 'B_upper', 0);
mustBeSize(data.c_upper, [n_x, 1], 'c_upper', 0);
mustBeSize(data.d_upper, [n_y, 1], 'd_upper', 0);
mustBeSize(data.rhs_upper, [n_G, 1], 'rhs_upper', 0);
% Lower-level arrays
mustBeSize(data.A_lower, [n_g, n_x], 'A_lower', 0);
mustBeSize(data.B_lower, [n_g, n_y], 'B_lower', 0);
mustBeSize(data.c_lower, [n_x, 1], 'c_lower', 0);
mustBeSize(data.d_lower, [n_y, 1], 'd_lower', 0);
mustBeSize(data.rhs_lower, [n_g, 1], 'rhs_lower', 0);
% Bounds (allow empty or correct length)
mustBeSize(data.x_lower_bound, [n_x, 1], 'x_lower_bound', 1);
mustBeSize(data.x_upper_bound, [n_x, 1], 'x_upper_bound', 1);
mustBeSize(data.y_lower_bound, [n_y, 1], 'y_lower_bound', 1);
mustBeSize(data.y_upper_bound, [n_y, 1], 'y_upper_bound', 1);
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=- %
%% Read JSON %%
% -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=- %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
function data = read_json(path)
% Read JSON file
txt = fileread(path);
data = jsondecode(txt);
% Store the path. This is useful diagnostics
data.path = char(path);
% ================== REQUIRED keys ==================
% Dimensions and objective coefficients are nessasary
keys = {
'n_x', 'n_y', 'c_upper', 'd_upper', 'c_lower', 'd_lower'
};
for i = 1:length(keys)
key = keys{i};
msg = sprintf('Missing required key "%s". See JSON file: %s.', key, path);
assert(isfield(data, key), msg)
end
% ================== OPTIONAL keys ==================
% Bounds are linear ineq constraints are optional
% They should be cast to the double data type
% Or default to the correct empty shape
keys = {
"A_upper"; "B_upper"; "c_upper"; "d_upper"; "rhs_upper"; "x_upper_bound"; "x_lower_bound"; ...
"A_lower"; "B_lower"; "c_lower"; "d_lower"; "rhs_lower"; "y_upper_bound"; "y_lower_bound";...
};
default_values = {
zeros(0, data.n_x); zeros(0, data.n_y); zeros(0, 1); zeros(0, 1); zeros(0, 1); zeros(0, 1); zeros(0, 1); ...
zeros(0, data.n_x); zeros(0, data.n_y); zeros(0, 1); zeros(0, 1); zeros(0, 1); zeros(0, 1); zeros(0, 1); ...
};
for i = 1:length(keys)
key = keys{i};
default_value = default_values{i};
if ~isfield(data, key) || isempty(data.(key))
data.(key) = default_value;
else
data.(key) = double(data.(key));
end
end
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=- %
%% Read AUX MPS Pair %%
% -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=- %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
function data = read_aux_mps(aux_path, mps_path)
% _READ_AUX_MPS
% Uses Gurobi MATLAB 'gurobi_read' for MPS,
% And parses the AUX file to separate UL/LL variables & constraints.
%
% NOTE: Requires Gurobi MATLAB interface on path.
arguments
aux_path string
mps_path string
end
% ---------- Read MSP ----------
model = gurobi_read(char(mps_path)); %#ok<NASGU> % returns struct
model = gurobi_read(char(mps_path)); % (MATLAB sometimes warns w/ first call)
% ---------- Read AUX ----------
aux = linear_bilevel.local_parse_aux_file(aux_path);
[var_names_lower, var_coeffs_lower] = linear_bilevel.local_parse_aux_variables(aux.VARS);
% ---------- Names ----------
var_names_all = string(model.varnames(:));
var_names_lower = string(var_names_lower);
var_names_upper = var_names_all(~ismember(var_names_all, var_names_lower));
constr_names_all = string(model.constrnames(:));
constr_names_lower = string(aux.CONSTRS(:));
constr_names_upper = constr_names_all(~ismember(constr_names_all, constr_names_lower));
% ---------- Dimentions ----------
n_y = length(var_names_lower);
n_x = length(var_names_upper);
n_constr_lower = numel(constr_names_lower);
n_constr_upper = numel(constr_names_upper);
% ---------- Indecies ----------
% Map names -> index in MPS columns
idx_lower = zeros(1, n_y);
idx_upper = zeros(1, n_y);
for i = 1:n_y
idx = find(var_names_all == var_names_lower(i), 1);
assert(~isempty(idx), "Lower-level variable '%s' not found in MPS.", var_names_lower(i));
idx_lower(i) = idx;
end
for i = 1:n_x
idx = find(var_names_all == var_names_upper(i), 1);
assert(~isempty(idx), "Upper-level variable '%s' not found in MPS.", var_names_upper(i));
idx_upper(i) = idx;
end
% Validation
assert(str2double(aux.NUMVARS) == n_y, "AUX NUMVARS=%d but @VARSBEGIN ... @VARSEND block has %d entries.", aux.NUMVARS, n_y);
% ---------- Allocate data ----------
data = struct( ...
'name', char(extractBefore(string(aux_path), ".aux")), ...
'path', char(aux_path), ...
'n_x', n_x, ...
'n_y', n_y, ...
'A_upper', zeros(n_constr_upper, n_x), ...
'B_upper', zeros(n_constr_upper, n_y), ...
'rhs_upper', zeros(n_constr_upper, 1), ...
'c_upper', zeros(n_x, 1), ...
'd_upper', zeros(n_y, 1), ...
'A_lower', zeros(n_constr_lower, n_x), ...
'B_lower', zeros(n_constr_lower, n_y), ...
'rhs_lower', zeros(n_constr_lower, 1), ...
'c_lower', zeros(n_x, 1), ...
'd_lower', var_coeffs_lower(:), ... % LL objective from AUX
'x_lower_bound', -inf(n_x, 1), ...
'x_upper_bound', inf(n_x, 1), ...
'y_lower_bound', -inf(n_y, 1), ...
'y_upper_bound', inf(n_y, 1) ...
);
% ---------- Constraints ----------
% Ax + By >= rhs
A = model.A; % sparse (m × n)
rhs = model.rhs(:); % (m × 1)
sense = string(model.sense(:)); % each of {'<','>','='} (or possibly 'L','G','E')
% Build name -> row index maps for upper/lower blocks
mapRowU = containers.Map(cellstr(constr_names_upper), num2cell(1:n_constr_upper));
mapRowL = containers.Map(cellstr(constr_names_lower), num2cell(1:n_constr_lower));
for i = 1:numel(constr_names_all)
cname = constr_names_all(i);
si = sense(i);
switch si
case "<"
sgn = -1;
case ">"
sgn = +1;
otherwise
error("Equality constraint '%s' detected. Only <= or >= are supported.", cname);
end
row_full = full(A(i, :));
rhs_i = sgn * rhs(i);
coeff_x = sgn * row_full(idx_upper);
coeff_y = sgn * row_full(idx_lower); % reordered to AUX (y) order
if isKey(mapRowU, cname)
r = mapRowU(cname);
data.A_upper(r, :) = coeff_x;
data.B_upper(r, :) = coeff_y;
data.rhs_upper(r) = rhs_i;
elseif isKey(mapRowL, cname)
r = mapRowL(cname);
data.A_lower(r, :) = coeff_x;
data.B_lower(r, :) = coeff_y;
data.rhs_lower(r) = rhs_i;
else
error("Constraint '%s' not found in either UL or LL sets.", cname);
end
end
% ---------- Lower-level objective ----------
% minimize c'x + d'y (from MPS)
obj = model.obj(:);
if isfield(model, 'modelsense')
if string(model.modelsense) == "min"
obj_sgn = +1;
elseif string(model.modelsense) == "max"
obj_sgn = -1;
else
error("Unknown model.modelsense: '%s'.", modelsense);
end
else
obj_sgn = +1;
end
data.c_upper = obj_sgn * obj(idx_upper);
data.d_upper = obj_sgn * obj(idx_lower);
% ---------- Bounds ----------
lb = model.lb(:);
ub = model.ub(:);
data.x_lower_bound = lb(idx_upper);
data.x_upper_bound = ub(idx_upper);
data.y_lower_bound = lb(idx_lower);
data.y_upper_bound = ub(idx_lower);
end
function aux = local_parse_aux_file(aux_path)
% Please see documentation at https://coin-or.github.io/MibS/input.html
%
% @VARSBEGIN ... @VARSEND (lines: "<name> <coeff>")
% @CONSTRSBEGIN ... @CONSTRSEND (lines: "<constraint_name>")
% @NUMVARS (single line integer on the next line)
raw = fileread(aux_path);
lines = regexp(raw, '\r\n|\n|\r', 'split');
lines = string(lines(:));
lines = strtrim(lines);
lines = lines(lines ~= "" & ~startsWith(lines, "#"));
aux = struct();
i = 1;
n = numel(lines);
% Process line by line
while i <= n
line = lines(i);
% Block sections.
% @VARSBEGIN or @CONSTRSBEGIN
tok = regexp(line, '^@([A-Z]+)\s*BEGIN$', 'tokens', 'once');
if ~isempty(tok)
section = string(tok{1});
i = i + 1;
block = strings(0,1);
% Read until end.
% @VARSEND or @CONSTRSEND
while i <= n && ~strcmpi(lines(i), "@" + section + "END")
block(end+1,1) = lines(i); %#ok<AGROW>
i = i + 1;
end
if i > n
error("AUX parse error: Missing '@%s END'.", section);
end
aux.(section) = block;
i = i + 1;
continue
end
% Single-line tokens.
% @NUMVARS, @NUMCONSTRS, @NAME, @MPS, @LP
tokens = ["@NUMVARS", "@NUMCONSTRS", "@NAME", "@MPS", "@LP"];
if ismember(line, tokens)
if i+1 > n
error("AUX parse error: Expected a value after '%s'.", line);
end
key = eraseBetween(line,1,1);
aux.(key) = lines(i+1);
i = i + 2;
continue
end
% Otherwise unrecognized line; skip/ignore (or error out)
% We ignore unknown directives to be permissive.
i = i + 1;
end
% Sanity defaults if missing blocks
assert(isfield(aux, 'VARS'), 'Whiles reading the AUX file, failed to find the token @VARSBEGIN ... @VARSEND')
assert(isfield(aux, 'CONSTRS'), 'Whiles reading the AUX file, failed to find the token @CONSTRSBEGIN ... @CONSTRSEND')
assert(isfield(aux, 'NUMVARS'), 'Whiles reading the AUX file, failed to find the token @NUMVARS')
end
function [names, coeffs] = local_parse_aux_variables(var_lines)
% LOCAL_PARSE_AUX_VARIABLES
% Input: string array like ["alpha 1.23", "beta 4.56", ...]
% Output: names (string array), coeffs (double column)
names = strings(0,1);
coeffs = zeros(0,1);
for k = 1:numel(var_lines)
parts = split(strtrim(var_lines(k)));
if numel(parts) ~= 2
error("Invalid VARS entry: '%s' (expected '<name> <coef>').", var_lines(k));
end
nm = parts(1);
val = str2double(parts(2));
if ~isfinite(val)
error("Invalid value (should be numeric) in VARS entry: '%s'.", var_lines(k));
end
names(end+1,1) = nm; %#ok<AGROW>
coeffs(end+1,1) = val; %#ok<AGROW>
end
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=- %
%% Dimension %%
% -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=- %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
function n = dimension(key, data)
% Key is the function/variable name
% Value is it's dimension
n = dictionary( ...
'x', data.n_x, ...
'y', data.n_y, ...
'F', 1, ...
'G', length(data.rhs_upper) + length(data.x_upper_bound) + length(data.x_lower_bound), ...
'H', 0, ...
'f', 1, ...
'g', length(data.rhs_lower) + length(data.y_upper_bound) + length(data.y_lower_bound), ...
'h', 0 ...
);
if isKey(n,key)
n = n(key);
end
end
end
endimport json
import os.path
import gurobipy as gp
from bolib3 import np
"""
Linear Bilevel Program (LBP)
============================
The LBP, sometimes called fully-linear, has both linear objective functions and linear constraints.
Each parameterization is defined by vectors of coefficients stored in bolib3/data/linear_bilevel.
Upper-level
===========
minimise c_upper ' x + d_upper ' y
subject to A_upper @ x + B_upper @ y >= rhs_upper
x >= x_lower_bound
x <= x_upper_bound
y solves lower-level
Lower-level
===========
minimise c_lower ' x + d_lower ' y
subject to A_lower @ x + B_lower @ y >= rhs_lower
y >= y_lower_bound
y <= y_upper_bound
Table of Dimensions
===================
| Array | Shape |
| --------- | ---------- |
| A_upper | (n_G, n_x) |
| B_upper | (n_G, n_y) |
| c_upper | (n_x,) |
| d_upper | (n_y,) |
| rhs_upper | (n_G,) |
| A_lower | (n_g, n_x) |
| B_lower | (n_g, n_y) |
| c_lower | (n_x,) |
| d_lower | (n_y,) |
| rhs_lower | (n_g,) |
were n_x and n_y are the number of upper-level and lower-level decision variables respectively.
and n_G and n_g are the number of upper-level and lower-level inequality constraints respectively.
Example Python use case
=======================
from bolib3 import linear_bilevel
data = linear_bilevel.read_data("bolib3/data/linear_bilevel/anandalingham_white_1990.json")
x = np.array([16])
y = np.array([11])
print(linear_bilevel.g(x, y, data)) # [28. 12. 0. 0. 12. 11.]
Accepted data formats
=====================
* JSON: A key-value format where the key is a variable name and the value is a vector of coefficients.
See the function _read_json() for more details.
* AUX + MPS: The auxiliary (AUX) file lists names of the lower-level variables and constraints.
The Mathematical Programming System (MPS) file defines all constraints in a name-based system.
See the function _read_aux_mps() for more details.
"""
# Properties
name: str = "linear_bilevel"
category: str = "archetype"
subcategory: str = ""
author: str = "Samuel Ward"
# Methods
def F(x, y, data):
"""
Upper-level objective function
minimise F(x,y) := c'x + d'y
"""
c_upper = data['c_upper']
d_upper = data['d_upper']
return np.dot(c_upper, x) + np.dot(d_upper, y)
def G(x, y, data):
"""
Upper-level inequality constraints
Ax + By >= rhs
x >= x_lower_bound
x <= x_upper_bound
"""
# Model parameters
A_upper = data['A_upper']
B_upper = data['B_upper']
rhs_upper = data['rhs_upper']
x_upper_bound = data['x_upper_bound']
x_lower_bound = data['x_lower_bound']
# Inequality constraints
ineq_list = [np.empty(0)]
if rhs_upper.size:
ineq_list.append(A_upper@x + B_upper@y - rhs_upper)
if x_lower_bound.size:
ineq_list.append(x - x_lower_bound)
if x_upper_bound.size:
ineq_list.append(x_upper_bound - x)
# Constraints G(x, y) >= 0
return np.concatenate(ineq_list)
def H(x, y, data=None):
"""
Upper-level equality constraints
"""
return np.empty(0)
def f(x, y, data):
"""
Lower-level objective function
minimise f(x,y) := c'x + d'y
"""
c_lower = data['c_lower']
d_lower = data['d_lower']
return np.dot(c_lower, x) + np.dot(d_lower, y)
def g(x, y, data):
"""
Lower-level inequality constraints
Ax + By >= rhs
y >= y_lower_bound
y <= y_upper_bound
"""
# Model parameters
A_lower = data['A_lower']
B_lower = data['B_lower']
rhs_lower = data['rhs_lower']
y_upper_bound = data['y_upper_bound']
y_lower_bound = data['y_lower_bound']
# Inequality constraints
ineq_constraints = [np.empty(0)]
if rhs_lower.size:
ineq_constraints.append(A_lower@x + B_lower@y - rhs_lower)
if y_lower_bound.size:
ineq_constraints.append(y - y_lower_bound)
if y_upper_bound.size:
ineq_constraints.append(y_upper_bound - y)
# Constraints g(x, y) >= 0
return np.concatenate(ineq_constraints)
def h(x, y, data=None):
"""
Lower-level equality constraints
"""
return np.empty(0)
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-
# Derivatives
# -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-
# =============================================================
def derivative(func, wrt=None):
"""
This is a mapping of strings to the public functions.
Rather than calling df_dx(x, y) you may use derivative("f", wrt="x")(x, y)
Rather than calling f(x, y) you may use derivative("f", wrt=None)(x, y)
"""
mapping = {
# Upper-level
"F": {None: F, "x": dF_dx, "y": dF_dy, "xx": d2F_dxx, "yy": d2F_dyy, "xy": d2F_dxy},
"G": {None: G, "x": dG_dx, "y": dG_dy},
"H": {None: H, "x": dH_dx, "y": dH_dy},
# Lower-level
"f": {None: f, "x": df_dx, "y": df_dy, "xx": d2f_dxx, "yy": d2f_dyy, "xy": d2f_dxy},
"g": {None: g, "x": dg_dx, "y": dg_dy},
"h": {None: h, "x": dh_dx, "y": dh_dy},
}
return mapping[func][wrt]
def dF_dx(x, y, data):
"""First-order derivative of the upper-level objective function F(x,y) with respect to x"""
return data['c_upper']
def dF_dy(x, y, data):
"""First-order derivative of the upper-level objective function F(x,y) with respect to y"""
return data['d_upper']
def d2F_dxx(x, y, data):
"""Second-order derivative (Hessian) of the upper-level objective function F(x,y) w.r.t. x twice"""
return np.zeros((data['n_x'], data['n_x']))
def d2F_dyy(x, y, data):
"""Second-order derivative (Hessian) of the upper-level objective function F(x,y) w.r.t. y twice"""
return np.zeros((data['n_y'], data['n_y']))
def d2F_dxy(x, y, data):
"""Second-order derivative (Hessian) of the upper-level objective function F(x,y) w.r.t. x and y"""
return np.zeros((data['n_x'], data['n_y']))
def dG_dx(x, y, data):
"""Jacobian matrix of the upper-level inequality constraints G(x,y) with respect to x"""
jacobian_list = [np.empty((0, data['n_x']))]
if data['A_upper'].size:
jacobian_list.append(data['A_upper'])
if data['x_lower_bound'].size:
jacobian_list.append(np.identity(data['n_x']))
if data['x_upper_bound'].size:
jacobian_list.append(-1*np.identity(data['n_x']))
return np.concatenate(jacobian_list)
def dG_dy(x, y, data):
"""Jacobian matrix of the upper-level inequality constraints G(x,y) with respect to y"""
jacobian_list = [np.empty((0, data['n_y']))]
if data['B_upper'].size:
jacobian_list.append(data['B_upper'])
if data['x_lower_bound'].size:
jacobian_list.append(np.zeros((data['n_x'], data['n_y'])))
if data['x_upper_bound'].size:
jacobian_list.append(np.zeros((data['n_x'], data['n_y'])))
return np.concatenate(jacobian_list)
def dH_dx(x, y, data=None):
"""Jacobian matrix of the upper-level equality constraints H(x,y) with respect to x"""
return np.empty(0)
def dH_dy(x, y, data=None):
"""Jacobian matrix of the upper-level equality constraints H(x,y) with respect to y"""
return np.empty(0)
def df_dx(x, y, data):
"""First-order derivative of the lower-level objective function f(x,y) with respect to x"""
return data['c_lower']
def df_dy(x, y, data):
"""First-order derivative of the lower-level objective function f(x,y) with respect to y"""
return data['d_lower']
def d2f_dxx(x, y, data):
"""Second-order derivative (Hessian) of the lower-level objective function f(x,y) with respect x twice"""
return np.zeros((data['n_x'], data['n_x']))
def d2f_dyy(x, y, data):
"""Second-order derivative (Hessian) of the lower-level objective function f(x,y) with respect to y twice"""
return np.zeros((data['n_y'], data['n_y']))
def d2f_dxy(x, y, data):
"""Second-order derivative (Hessian) of the lower-level objective function f(x,y) with respect to x and y"""
return np.zeros((data['n_x'], data['n_y']))
def dg_dx(x, y, data):
"""Jacobian matrix of the lower-level inequality constraints g(x,y) with respect to x"""
jacobian_list = [np.empty((0, data['n_x']))]
if data['A_lower'].size:
jacobian_list.append(data['A_lower'])
if data['y_lower_bound'].size:
jacobian_list.append(np.zeros((data['n_y'], data['n_x'])))
if data['y_upper_bound'].size:
jacobian_list.append(np.zeros((data['n_y'], data['n_x'])))
return np.concatenate(jacobian_list)
def dg_dy(x, y, data):
"""Jacobian matrix of the lower-level inequality constraints g(x,y) with respect to y"""
jacobian_list = [np.empty((0, data['n_y']))]
if data['B_lower'].size:
jacobian_list.append(data['B_lower'])
if data['y_lower_bound'].size:
jacobian_list.append(np.identity(data['n_y']))
if data['y_upper_bound'].size:
jacobian_list.append(-1*np.identity(data['n_y']))
return np.concatenate(jacobian_list)
def dh_dx(x, y, data=None):
"""Jacobian matrix of the lower-level equality constraints h(x,y) with respect to x"""
return np.empty(0)
def dh_dy(x, y, data=None):
"""Jacobian matrix of the lower-level equality constraints h(x,y) with respect to y"""
return np.empty(0)
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-
# Read Data
# -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-
# =============================================================
json_datasets: list = sorted([
'anandalingham_white_1990.json', 'bard_1984a.json', 'bard_1984b.json', 'bard_1991_ex2.json',
'bard_falk_1982_ex2.json', 'ben_ayed_blair_1990a.json', 'ben_ayed_blair_1990b.json', 'bialas_karwan_1984a.json',
'bialas_karwan_1984b.json', 'candler_townsley_1982.json', 'clark_westerberg_1988.json',
'clark_westerberg_1990b.json', 'glackin_et_al_2009.json', 'haurie_savard_white_1990.json',
'hu_huang_zhang_2009.json', 'hypercube_10.json', 'hypercube_100.json', 'hypercube_3.json', 'hypercube_50.json',
'lan_wen_shih_lee_2007.json', 'liu_hart_1994.json', 'mersha_dempe_2006_ex1.json', 'mersha_dempe_2006_ex2.json',
'tuy_et_al_1993.json', 'tuy_et_al_1994.json', 'tuy_et_al_2007_ex3.json', 'visweswaran_et_al_1996.json',
'wang_jiao_li_2005.json'
])
aux_datasets: list = sorted([
'general20-20-10-20-20-1.aux', 'BCPIns_8_7_1.txt_trad.txt_K5.aux'
])
datasets: list = json_datasets + aux_datasets
paths: list = (
[os.path.join('bolib3', 'data', 'linear_bilevel', dataset) for dataset in json_datasets] +
[os.path.join('bolib3', 'data', 'linear_bobilib_aux_mps', dataset) for dataset in aux_datasets]
)
def read_data(path: str = datasets[0]) -> dict:
"""A path to either a .json or .aux file"""
# JSON
if path.endswith('.json'):
data = _read_json(path)
# AUX
elif path.endswith('.aux'):
data_name = os.path.basename(path).split('.aux')[0]
directory = os.path.dirname(path)
# Find the corresponding MPS file
if os.path.exists(os.path.join(directory, data_name + '.mps')):
mps_path = os.path.join(directory, data_name + '.mps')
elif os.path.exists(os.path.join(directory, data_name + '.mps.gz')):
mps_path = os.path.join(directory, data_name + '.mps.gz')
else:
raise FileNotFoundError(f"Failed to find {data_name}.mps or {data_name}.mps.gz file in {directory}")
data = _read_aux_mps(aux_path=path, mps_path=mps_path, data_name=data_name)
# Error
else:
raise ValueError(f"Dataset for linear_bilevel should be either .json or .aux not '{path}'")
# Validate the dimensions of all arrays
_validate_data_dimensions(data)
return data
def _validate_data_dimensions(data: dict):
"""
See "Table of Dimensions" in the docstring at the top of this file.
"""
n_x, n_y, n_G, n_g = data['n_x'], data['n_y'], data['rhs_upper'].shape[0], data['rhs_lower'].shape[0]
identifier = "; ".join((data.get('name', ''), data.get('path', '')))
for key, expected_shape in {
'A_upper': (n_G, n_x),
'B_upper': (n_G, n_y),
'c_upper': (n_x,),
'd_upper': (n_y,),
'rhs_upper': (n_G,),
'A_lower': (n_g, n_x),
'B_lower': (n_g, n_y),
'c_lower': (n_x,),
'd_lower': (n_y,),
'rhs_lower': (n_g,),
}.items():
assert data[key].shape == expected_shape, \
f"Dimension of {key} is {data[key].shape} but was expect to be {expected_shape}. See: '{identifier}'."
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-
# Read JSON
# -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-
# =============================================================
def _read_json(path: str) -> dict:
"""
Must give a path to a JSON file with all the keys listed in required_data_fields.
For example, it should begin something like this:
{
"name": "author_1990",
"n_x": 2,
"n_y": 3,
"A_upper": [
[4.0, 5.0]
],
"""
# Load the JSON file
with open(path, 'r') as file:
data = json.load(file)
# Store the path. This is useful diagnostics
data['path'] = path
# ================== REQUIRED keys ==================
# Dimensions and objective coefficients are nessasary
# They must be cast to the correct type
for key, expected_type in {
'n_x': int, 'c_upper': np.array, 'd_upper': np.array,
'n_y': int, 'c_lower': np.array, 'd_lower': np.array,
}.items():
assert key in data, f"Missing required key \"{key}\". See JSON file: {path}."
try:
data[key] = expected_type(data[key])
except TypeError:
raise (f"Value of \"{key}\": {data[key]} should be type \'{expected_type.__name__}\'. "
f"See in JSON file: {path}.")
# ================== OPTIONAL keys ==================
# Bounds and linear ineq constraints are optional
# They should be cast to numpy arrays
# Or default to the correct empty shape
for key, default_shape in {
"x_upper_bound": (0,), "x_lower_bound": (0,), "y_upper_bound": (0,), "y_lower_bound": (0,),
'A_upper': (0, data['n_x']), 'B_upper': (0, data['n_y']), "rhs_upper": (0,),
'A_lower': (0, data['n_x']), 'B_lower': (0, data['n_y']), "rhs_lower": (0,),
}.items():
if (key in data) and data[key]:
data[key] = np.array(data[key])
else:
data[key] = np.empty(default_shape)
return data
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-
# Read AUX MPS pair
# -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-
# =============================================================
def _read_aux_mps(aux_path: str, mps_path: str, data_name: str) -> dict:
"""
Please see documentation at https://coin-or.github.io/MibS/input.html
:param aux_path: An auxiliary (aux) file that contains information necessary to
identify the upper and lower-level constraints and variables.
:param mps_path: An instance file in MPS, LP, or GMPL/AMPL format that contains a
description of the upper level objective, the variables (upper and
lower level), and the constraints (upper and lower level).
:return: dictionary of numpy arrays
"""
# Read MPS file
gurobi_model = gp.read(mps_path)
# Read AUX file
with open(aux_path, 'r') as file:
lines = [line.strip() for line in file if line.strip() and not line.strip().startswith("#")]
# //==================================================\\
# || Process AUX lines ||
# \\==================================================//
aux_data = {}
current_section = None
collecting = False
for line in lines:
if not line or line.startswith("#"):
continue
if line.startswith("@"):
if line in ["@NUMVARS", "@NUMCONSTRS", "@NAME", "@MPS", "@LP"]:
current_section = line[1:]
collecting = False
elif line.endswith("BEGIN"):
current_section = line[1:-5] # strip @ and BEGIN
aux_data[current_section] = []
collecting = True
elif line.endswith("END"):
collecting = False
current_section = None
continue
if current_section:
if collecting:
aux_data[current_section].append(line)
else:
# Single-line entries
aux_data[current_section] = line
# Assert that these keys must exist
for key, tocken_form in {
"VARS": "@VARSBEGIN ... @VARSEND",
"CONSTRS": "@CONSTRSBEGIN ... @CONSTRSEND",
"NUMVARS": "@NUMVARS"
}.items():
message = (f"Whiles reading the AUX file, failed to find the token '{tocken_form}'. "
f"\n See the documentation: https://coin-or.github.io/MibS/input.html"
f"\n Check the AUX file: '{aux_path}'.")
assert key in aux_data, message
# Lists of upper and lower variables and constraints
var_names_lower, var_coeffs_lower = _parse_aux_variables(aux_data.get("VARS", []))
var_names_upper = [v.VarName for v in gurobi_model.getVars() if v.VarName not in var_names_lower]
constr_names_lower = aux_data['CONSTRS']
constr_names_upper = [c.ConstrName for c in gurobi_model.getConstrs() if c.ConstrName not in aux_data['CONSTRS']]
# Dimensions
n_y = int(aux_data['NUMVARS'])
n_x = (len(gurobi_model.getVars()) - n_y)
n_constr_lower = len(aux_data['CONSTRS'])
n_constr_upper = len(gurobi_model.getConstrs()) - n_constr_lower
# Construct the data dictionary
data = {
"name": data_name,
"n_x": n_x,
"n_y": n_y,
"A_upper": np.zeros((n_constr_upper, n_x)),
"B_upper": np.zeros((n_constr_upper, n_y)),
"rhs_upper": np.zeros(n_constr_upper),
"c_upper": np.zeros(n_x),
"d_upper": np.zeros(n_y),
"A_lower": np.zeros((n_constr_lower, n_x)),
"B_lower": np.zeros((n_constr_lower, n_y)),
"rhs_lower": np.zeros(n_constr_lower),
"c_lower": np.zeros(n_x),
"d_lower": np.array(var_coeffs_lower),
"x_lower_bound": -np.inf*np.ones(n_x),
"x_upper_bound": np.inf*np.ones(n_x),
"y_lower_bound": -np.inf*np.ones(n_y),
"y_upper_bound": np.inf*np.ones(n_y),
}
# //==================================================\\
# || Constraints ||
# \\==================================================//
# Coefficients of the linear inequality constraints
# Ax + By >= rhs
# Iterate over all the MPS constraints
for gurobi_constr in gurobi_model.getConstrs():
constr_name = gurobi_constr.ConstrName
constr_expr = gurobi_model.getRow(gurobi_constr)
# Multiply by -1 to flip inequality
match gurobi_constr.Sense:
case '<' | '<=' | r'\leq':
constr_sign = -1
case '>' | '>=' | r'\geq':
constr_sign = +1
case _:
raise ValueError(
f"For an inequality constraint, gurobi_constr.Sense should be '<' or '>' not {gurobi_constr.Sense}")
# Right hand side
if constr_name in constr_names_upper:
level = 'upper'
row = constr_names_upper.index(constr_name)
elif constr_name in constr_names_lower:
level = 'lower'
row = constr_names_lower.index(constr_name)
else:
raise ValueError(f"Failed to find constraint '{constr_name}'.")
data[f'rhs_{level}'][row] = constr_sign*gurobi_constr.RHS
# Iterate over the terms of the constraint
for i in range(constr_expr.size()):
var_name = constr_expr.getVar(i).VarName
coeff = constr_expr.getCoeff(i)
if var_name in var_names_upper:
matrix = 'A'
col = var_names_upper.index(var_name)
elif var_name in var_names_lower:
matrix = 'B'
col = var_names_lower.index(var_name)
else:
raise ValueError(f"Failed to find variable '{var_name}'.")
data[f"{matrix}_{level}"][row, col] = constr_sign*coeff
# //==================================================\\
# || Objective ||
# \\==================================================//
# Objective coefficients
# minimise c'x + d'y
# Extract from gurobi
obj_expr = gurobi_model.getObjective()
obj_sign = gurobi_model.ModelSense # 1 for minimization, -1 for maximization
for i in range(obj_expr.size()):
var_name = obj_expr.getVar(i).VarName
var_coeff = obj_expr.getCoeff(i)
if var_name in var_names_upper:
idx = var_names_upper.index(var_name)
data['c_upper'][idx] = obj_sign*var_coeff
elif var_name in var_names_lower: # upper level var
idx = var_names_lower.index(var_name)
data['d_upper'][idx] = obj_sign*var_coeff
else:
raise ValueError(f"Variable '{var_name}' not found in either LL or UL variable lists.")
# //==================================================\\
# || Bounds ||
# \\==================================================//
for var in gurobi_model.getVars():
var_name, lb, up = var.VarName, var.LB, var.UB
if var_name in var_names_upper:
idx = var_names_upper.index(var_name)
data['x_lower_bound'][idx] = lb
data['x_upper_bound'][idx] = up
elif var_name in var_names_lower:
idx = var_names_lower.index(var_name)
data['y_lower_bound'][idx] = lb
data['y_upper_bound'][idx] = up
else:
raise ValueError(f"Variable '{var_name}' not found in UL or LL variable lists.")
# Return the dictionary
return data
def _parse_aux_variables(
List_of_var_strings: list[str]
) -> tuple[list[str], list[float]]:
"""
Converts a list of strings of the form: ["{name} {coefficient}", ...]
into two separate lists, one of names and one of coefficients
For example:
List_of_var_strings := ["alpha 1.23", "beta 4.56", "gamma 7.89", ...]
list_of_name := ["alpha", "beta", "gamma" ...]
list_of_coefficients := [1.23, 4.56, 7.89, ...]
"""
list_of_name: list[str] = []
list_of_coefficients: list[float] = []
for var in List_of_var_strings:
parts = var.split()
if len(parts) != 2:
raise ValueError(f"Invalid VARS entry: {var}")
try:
var_name = str(parts[0])
value = float(parts[1])
except ValueError:
raise ValueError(f"Invalid numeric value in VARS entry: {var}")
list_of_name.append(var_name)
list_of_coefficients.append(value)
return list_of_name, list_of_coefficients
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-
# Dimensions
# -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-
# =============================================================
def dimension(key='', data=None):
"""
If the argument 'key' is not specified, then:
- a dictionary mapping variable/function names (str) to the corresponding dimension (int) is returned.
If the first argument 'key' is specified, then:
- a single integer representing the dimension of the variable/function with the name {key} is returned.
"""
n = {
"x": data['n_x'],
"y": data['n_y'],
"F": 1,
"G": len(data['rhs_upper']) + len(data['x_upper_bound']) + len(data['x_lower_bound']),
"H": 0,
"f": 1,
"g": len(data['rhs_lower']) + len(data['y_upper_bound']) + len(data['y_lower_bound']),
"h": 0,
}
if key in n:
return n[key]
return n