Added 15/09/2025
Archetype
polynomial_bilevel
Datasets
Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 3,
"H": 0,
"f": 1,
"g": 3,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 4,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Description
Dempe, Stephan and Dutta, Joydeep (2012) (see page 42, Example 2.4).
Is bilevel programming a special case of a mathematical program with complementarity constraints?.
https://doi.org/10.1007/s10107-010-0342-1Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 0,
"H": 0,
"f": 1,
"g": 1,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 6,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 3,
"H": 0,
"f": 1,
"g": 3,
"h": 0
}Solution
{
"optimality": "unknown"
}Description
Ishizuka, Yo and Aiyoshi, Eitaro (1992) (see page 82, Example 1).
Double penalty method for bilevel optimization problems.
https://doi.org/10.1007/BF02098173Dimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 5,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Description
Lampariello, Lorenzo and Sagratella, Simone (2020) (see page 281, Example 2).
Numerically tractable optimistic bilevel problems.
https://doi.org/10.1007/s10589-020-00178-yDimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.10).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [1],
"y": [0.5],
"F": 0.5,
"G": [0.9,0],
"H": [],
"f": -1,
"g": [1.5,0.5],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.11).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [0],
"y": [-0.8],
"F": -0.8,
"G": [1,1],
"H": [],
"f": 0,
"g": [0,1.8],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.12).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [0],
"y": [0],
"F": 0,
"G": [1,1],
"H": [],
"f": 0,
"g": [1,1],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.13).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [0],
"y": [1],
"F": -1,
"G": [1,1],
"H": [],
"f": 0,
"g": [2,0],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.14).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [0.25],
"y": [0.5],
"F": 0.25,
"G": [1.25,0.75],
"H": [],
"f": -0.08333333333333334,
"g": [1.5,0.5],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.15).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [-1],
"y": [1],
"F": 0,
"G": [0,2],
"H": [],
"f": -0.8333333333333333,
"g": [2,0],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.16).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [-0.5],
"y": [-1],
"F": -2,
"G": [0.5,1.5],
"H": [],
"f": 0,
"g": [0,2],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.17).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [-0.25],
"y": [0.5],
"F": 0.1875,
"G": [0.75,1.25],
"H": [],
"f": -0.015625,
"g": [1.5,0.5],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.18).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [0.5],
"y": [0],
"F": -0.25,
"G": [1.5,0.5],
"H": [],
"f": 0,
"g": [1,1],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.19).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [0.18858048469644506],
"y": [-0.43425854591066493],
"F": 0.44665610118748084,
"G": [1.1885804846964452,0.8114195153035549],
"H": [],
"f": -0.017781299604173073,
"g": [0.5657414540893351,1.434258545910665],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.20).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [0.5],
"y": [0.5],
"F": 0.3125,
"G": [1.5,0.5],
"H": [],
"f": -0.08333333333333334,
"g": [1.5,0.5],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.21).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [-0.5545],
"y": [0.4554],
"F": 1.54025941,
"G": [0.4455,1.5545],
"H": [],
"f": -0.0474164946863091,
"g": [1.4554,0.5446],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.22).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 3,
"h": 0
}Solution
{
"optimality": "global",
"x": [-0.5545],
"y": [0.4554],
"F": 1.54025941,
"G": [0.4455,1.5545],
"H": [],
"f": -0.0474164946863091,
"g": [0.1943144575,1.4554,0.5446],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.23).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 3,
"H": 0,
"f": 1,
"g": 3,
"h": 0
}Solution
{
"optimality": "global",
"x": [-0.5545],
"y": [0.4554],
"F": 0.30747025,
"G": [2.77713225,0.4455,1.5545],
"H": [],
"f": 0.4554,
"g": [0.21869186922,1.4554,0.5446],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.24).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.25).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 2,
"y": 3,
"F": 1,
"G": 6,
"H": 0,
"f": 1,
"g": 9,
"h": 0
}Solution
{
"optimality": "best_known",
"x": [-1,-1],
"y": [-1,1,1],
"F": -1,
"G": [1.1,1,0,0,2,2],
"H": [],
"f": -2,
"g": [0,0,0.9,0,2,2,2,0,0],
"h": []
}Description
Mitsos, Alexander and Barton, Paul I (2006) (see Example 3.26).
A test set for bilevel programs.
https://www.researchgate.net/publication/228455291Dimension
{
"x": 2,
"y": 3,
"F": 1,
"G": 7,
"H": 0,
"f": 1,
"g": 6,
"h": 0
}Solution
{
"optimality": "global",
"x": [-1,-1],
"y": [1,1,-0.7071067811865476],
"F": -2.353553390593274,
"G": [0.9,1,-1.1102230246251565e-16,0,0,2,2],
"H": [],
"f": -2,
"g": [2,2,0.2928932188134524,0,0,1.7071067811865475],
"h": []
}Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 3,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Description
Nie, Jiawang and Wang, Li and Ye, Jane J. (2017) (see page 1740, Example 3.4).
Bilevel Polynomial Programs and Semidefinite Relaxation Methods.
https://doi.org/10.1137/15M1052172Dimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Description
Nie, Jiawang and Wang, Li and Ye, Jane J. (2017) (see page 1747, Example 5.2).
Bilevel Polynomial Programs and Semidefinite Relaxation Methods.
https://doi.org/10.1137/15M1052172Dimension
{
"x": 2,
"y": 3,
"F": 1,
"G": 5,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Description
Nie, Jiawang and Wang, Li and Ye, Jane J. (2017) (see page 1748, Example 5.4).
Bilevel Polynomial Programs and Semidefinite Relaxation Methods.
https://doi.org/10.1137/15M1052172Dimension
{
"x": 4,
"y": 4,
"F": 1,
"G": 3,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Description
Nie, Jiawang and Wang, Li and Ye, Jane J. (2017) (see page 1750, Example 5.7).
Bilevel Polynomial Programs and Semidefinite Relaxation Methods.
https://doi.org/10.1137/15M1052172Dimension
{
"x": 2,
"y": 3,
"F": 1,
"G": 5,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Description
Nie, Jiawang and Wang, Li and Ye, Jane J. (2017) (see page 1751, Example 5.8).
Bilevel Polynomial Programs and Semidefinite Relaxation Methods.
https://doi.org/10.1137/15M1052172Dimension
{
"x": 4,
"y": 4,
"F": 1,
"G": 3,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Description
Nie, Jiawang and Wang, Li and Ye, Jane J. (2017) (see page 1753, Example 6.1).
Bilevel Polynomial Programs and Semidefinite Relaxation Methods.
https://doi.org/10.1137/15M1052172Dimension
{
"x": 2,
"y": 2,
"F": 1,
"G": 5,
"H": 0,
"f": 1,
"g": 1,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 4,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 4,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 2,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 4,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 4,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 4,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 1,
"h": 0
}Solution
{
"optimality": "unknown"
}Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 2,
"H": 0,
"f": 1,
"g": 1,
"h": 0
}Solution
{
"optimality": "unknown"
}\subsection{polynomial\_bilevel}
\label{subsec:polynomial_bilevel}
% Polynomial in x and y
\newcommand{\poly}[1]{\mathcal{P}_{#1}}
Let $\x\in\R^{n_x},\y\in\R^{n_y}$ be the lower and upper level decision variables, respectively.
Given a data parametrisation $m\in\N$, the number of terms; $c\in\R^m$ a vector of coefficients; $p\in\Z^{n_x,m}$ a matrix of exponents for $x$; and $p\in\Z^{m,n_y}$ a matrix of exponents for $y$ then we define the polynomial as:
\begin{equation*}
\poly{c,p,q}(x,y)\coloneqq \sum_{K=1}^{m} c_k \prod^{n_x}_{i=1} x_i^{p_{k,i}}\prod^{n_y}_{j=1}y^{q_{k,j}}.
\end{equation*}
% Polynomial Bilevel Program (PBP)
We make no restriction on the degree or number of terms in these polynomials.
Given a sequence of polynomials $\poly{\upper}^0, \poly{\upper}^k, \poly{\low}^0, \poly{\low}^k$, the The \emph{Polynomial Bilevel Program} \eqref{eq:PBP} is defined as:
\begin{flalign*}
\label{eq:PBP}
\tag{PBP}
\minimise_{\x\in\R^{n_x}, \y\in\R^{n_y}} \quad
& \poly{\upper}^0(x,y)\\
\subjectto \quad
& \poly{\upper}^k(x,y)\geq 0 \qquad \text{ for } k=1,2,\dots,\\
& {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}
\poly{\low}^0(x,y)\\
\st[1]\quad \poly{\low}^k(x,y)\geq 0, \qquad \text{ for } k=1,2,\dots,\\
\st[2]\quad {A_\low} \x + {B_\low} \y \geq {r_\low},\\
\st[3]\quad \ylb\leq \y \leq\yub.
\end{array}
\right.
\end{flalign*}
% Demonstrate JSON
The~\eqref{eq:PBP} is represented neatly by a JSON file.
Below, we demonstrate how the lower-level objective function of Mitsos and Barton's test example~\cite[Example 3.10]{Mitsos2006} is represented in JSON by the parameters $(m, c, p, q)$ of its polynomial.
\begin{align*}
\begin{array}{llllll}
f(x,y)=&16xy^4&+2xy^3&-8xy^2&-1.5xy&+0.5x.
\end{array}
\end{align*}
\begin{lstlisting}[style=custom,language=JSON,label={lst:polynomial_json}]
"f": {
"m": 5,
"c": [ 16.0, 2.0, -8.0, -1.5, 0.5 ],
"p": [[ 1, 1, 1, 1, 1 ]],
"q": [[ 4, 3, 2, 1, 0 ]]
}
\end{lstlisting}
% Remarks
This allows for easy derivative calculations using the formal derivative.
Furthermore, the work of Nie et al.~\cite{Nie2017,Nie2021} suggests a bespoke numerical method for solving the polynomial bilevel program.classdef polynomial_bilevel
% Polynomial Bilevel Program (PBP)
% ================================
% The objective function and constraints are polynomials.
% Each polynomial is parameterised by m, c, p, q as defined below.
% These parameters can be loaded from the JSON files in bolib3/data/polynomial_bilevel.
%
%
% Polynomial expression
% =====================
% m: Number of terms.
% c: Vector (m,) of coefficients.
% p: Matrix (m, n_x) of exponents for x.
% q: Matrix (m, n_y) of exponents for y.
%
% In LaTeX syntax:
% Polynomial(x, y) := sum_{t=1,...,m}(
% c_t
% * product_{i=1,...,n_x}( x_i^{p_{t,i}} )
% * product_{i=1,...,n_y}( y_i^{q_{t,i}} )
% )
%
% In MATLAB syntax:
% Polynomial(x, y) := sum(c .* prod(x'.^p, 2) .* prod(y'.^q, 2))
properties(Constant)
name = 'polynomial_bilevel';
category = 'archetype';
subcategory = '';
author = 'Samuel Ward';
datasets = {
'dempe_dutta_2012_ex24.json'; 'colson_2002_bipa1.json'; 'colson_2002_bipa2.json'; 'colson_2002_bipa4.json';
'floudas_zlobec_1998.json'; 'gumus_floudas_2001_ex1.json'; 'ishizuka_aiyoshi_1992a.json';
'kleniati_adjiman_2014_ex3.json'; 'lampariello_sagratella_2017_ex23.json';
'mitsos_barton_2006_ex39.json'; 'mitsos_barton_2006_ex310.json'; 'mitsos_barton_2006_ex311.json';
'mitsos_barton_2006_ex312.json'; 'mitsos_barton_2006_ex313.json'; 'mitsos_barton_2006_ex314.json';
'mitsos_barton_2006_ex315.json'; 'mitsos_barton_2006_ex316.json'; 'mitsos_barton_2006_ex317.json';
'mitsos_barton_2006_ex318.json'; 'mitsos_barton_2006_ex319.json'; 'mitsos_barton_2006_ex320.json';
'mitsos_barton_2006_ex321.json'; 'mitsos_barton_2006_ex322.json'; 'mitsos_barton_2006_ex323.json';
'mitsos_barton_2006_ex324.json'; 'mitsos_barton_2006_ex325.json'; 'mitsos_barton_2006_ex326.json';
'nie_wang_ye_2017_ex34.json'; 'nie_wang_ye_2017_ex52.json'; 'nie_wang_ye_2017_ex54.json';
'nie_wang_ye_2017_ex57.json'; 'nie_wang_ye_2017_ex58.json'; 'nie_wang_ye_2017_ex61.json';
'outrata_1993_ex31.json'; 'outrata_1993_ex32.json'; 'outrata_1994_ex31.json';
'paulavicius_adjiman_2017a.json'; 'paulavicius_adjiman_2017b.json'; 'shimizu_etal_1997b.json';
'vogel_2012.json'; 'ye_zhu_2010_ex42.json'; 'ye_zhu_2010_ex43.json';
};
paths = fullfile('bolib3', 'data', 'polynomial_bilevel', polynomial_bilevel.datasets);
end
methods(Static)
function val = polynomial(x, y, coeff_expon)
% Polynomial in
% x: vector (n_x, 1)
% y: vector (n_x, 1)
% The coeff_expon argument should be a stuct with fields
% c: vector (m, 1) of coefficients.
% p: matrix (m, n_x) of exponents for x.
% q: matrix (m, n_y) of exponents for y.
val = sum(coeff_expon.c .* ...
prod(power(x(:).', coeff_expon.p), 2) .* ...
prod(power(y(:).', coeff_expon.q), 2) ...
);
end
function val = F(x, y, data)
% Upper-level objective function
val = polynomial_bilevel.polynomial(x, y, data.F);
end
function ineq = G(x, y, data)
% Upper-level inequality constraints
ineq = [];
if ~isempty(data.ineq_upper)
ineq = zeros(length(data.ineq_upper), 1);
for i = 1:length(data.ineq_upper)
ineq(i) = polynomial_bilevel.polynomial(x, y, data.ineq_upper(i));
end
end
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 val = f(x, y, data)
% Lower-level objective function
val = polynomial_bilevel.polynomial(x, y, data.f);
end
function val = g(x, y, data)
% Lower-level inequality constraints
val = [];
if ~isempty(data.ineq_lower)
val = zeros(length(data.ineq_lower), 1);
for i = 1:length(data.ineq_lower)
val(i) = polynomial_bilevel.polynomial( x, y, data.ineq_lower(i));
end
end
if ~isempty(data.rhs_lower)
val = [val; data.A_lower * x + data.B_lower * y - data.rhs_lower];
end
if ~isempty(data.y_lower_bound)
val = [val; y - data.y_lower_bound];
end
if ~isempty(data.y_upper_bound)
val = [val; 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(x, y, data)
% ∂F/∂x
% Dimension (n_x × 1)
dx = polynomial_bilevel.polynomial_derivative(x, y, data.F, 'x');
end
function dy = dF_dy(x, y, data)
% ∂F/∂y
% Dimension (n_y × 1)
dy = polynomial_bilevel.polynomial_derivative(x, y, data.F, 'y');
end
function Hxx = d2F_dxx(x, y, data)
% ∂²F/∂x²
% Dimension (n_x × n_x)
Hxx = polynomial_bilevel.polynomial_derivative(x, y, data.F, 'xx');
end
function Hyy = d2F_dyy(x, y, data)
% ∂²F/∂y²
% Dimension (n_y × n_y)
Hyy = polynomial_bilevel.polynomial_derivative(x, y, data.F, 'yy');
end
function Hxy = d2F_dxy(x, y, data)
% ∂²F/∂x∂y
% Dimension (n_x × n_y)
Hxy = polynomial_bilevel.polynomial_derivative(x, y, data.F, 'xy');
end
function Jx = dG_dx(x, y, data)
% Jacobian of G wrt x
% Dimension (|G| × n_x)
Jx = zeros(0, data.n_x);
if ~isempty(data.ineq_upper)
Jx = zeros(length(data.ineq_upper), data.n_x);
for i = 1:length(data.ineq_upper)
Jx(i,:) = polynomial_bilevel.polynomial_derivative( ...
x, y, data.ineq_upper(i), 'x' ...
);
end
end
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(x, y, data)
% Jacobian of G wrt y:
% Dimension (|G| × n_y)
Jy = zeros(0, data.n_y);
if ~isempty(data.ineq_upper)
Jy = zeros(length(data.ineq_upper), data.n_y);
for i = 1:length(data.ineq_upper)
Jy(i,:) = polynomial_bilevel.polynomial_derivative( ...
x, y, data.ineq_upper(i), 'y' ...
);
end
end
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
% Dimension (0 × n_x)
Jx = zeros(0, data.n_x);
end
function Jy = dH_dy(~, ~, data)
% Jacobian of H wrt y
% Dimension (0 × n_y)
Jy = zeros(0, data.n_y);
end
function dx = df_dx(x, y, data)
% ∂f/∂x
% Dimension (n_x × 1)
dx = polynomial_bilevel.polynomial_derivative(x, y, data.f, 'x');
end
function dy = df_dy(x, y, data)
% ∂f/∂y
% Dimension (n_y × 1)
dy = polynomial_bilevel.polynomial_derivative(x, y, data.f, 'y');
end
function Hxx = d2f_dxx(x, y, data)
% ∂²f/∂x²
% Dimension (n_x × n_x)
Hxx = polynomial_bilevel.polynomial_derivative(x, y, data.f, 'xx');
end
function Hyy = d2f_dyy(x, y, data)
% ∂²f/∂y²
% Dimension (n_y × n_y)
Hyy = polynomial_bilevel.polynomial_derivative(x, y, data.f, 'yy');
end
function Hxy = d2f_dxy(x, y, data)
% ∂²f/∂x∂y
% Dimension (n_x × n_y)
Hxy = polynomial_bilevel.polynomial_derivative(x, y, data.f, 'xy');
end
function Jx = dg_dx(x, y, data)
% Jacobian of g wrt x
% Dimension (|g| × n_x)
Jx = zeros(0, data.n_x);
if ~isempty(data.ineq_lower)
Jx = zeros(length(data.ineq_lower), data.n_x);
for i = 1:length(data.ineq_lower)
Jx(i,:) = polynomial_bilevel.polynomial_derivative( ...
x, y, data.ineq_lower(i), 'x' ...
);
end
end
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(x, y, data)
% Jacobian of g wrt y
% Dimension (|g| × n_y)
Jy = zeros(0, data.n_y);
if ~isempty(data.ineq_lower)
Jy = zeros(length(data.ineq_lower), data.n_y);
for i = 1:length(data.ineq_lower)
Jy(i,:) = polynomial_bilevel.polynomial_derivative( ...
x, y, data.ineq_lower(i), 'y' ...
);
end
end
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
% Dimension (0 × n_x)
Jx = zeros(0, data.n_x);
end
function Jy = dh_dy(~, ~, data)
% Jacobian of h wrt y
% Dimension (0 × n_y)
Jy = zeros(0, data.n_y);
end
function out = polynomial_derivative(x, y, coeff_expon, with_respect_to)
% Compute first/second derivatives of a mixed polynomial in x and y.
%
% Inputs
% ------
% x : vector (n_x, 1)
% y : vector (n_y, 1)
% coeff_expon : struct with fields
% .c : (m, 1) vector of coefficients
% .p : (m, n_x) exponents for x
% .q : (m, n_y) exponents for y
% with_respect_to : char, one of 'x', 'y', 'xx', 'xy', 'yy'
% 'x' -> gradient wrt x (n_x, 1)
% 'y' -> gradient wrt y (n_y, 1)
% 'xx' -> Hessian wrt x (n_x, n_x)
% 'xy' -> cross Hessian (n_x, n_y)
% 'yy' -> Hessian wrt y (n_y, n_y)
%
% Output
% ------
% out : derivative (vector or matrix) as described above.
arguments
x double {mustBeVector}
y double {mustBeVector}
coeff_expon struct
with_respect_to char {mustBeMember(with_respect_to, {'x','y','xx','xy','yy'})}
end
% Normalize shapes
x_row = x(:).'; % (1, n_x)
y_row = y(:).'; % (1, n_y)
% Extract coefficients and exponents
c = coeff_expon.c(:); % (m x 1)
p = coeff_expon.p; % (m x n_x)
q = coeff_expon.q; % (m x n_y)
[~, n_x] = size(p); assert(n_x ~= 0);
[~, n_y] = size(q); assert(n_y ~= 0);
% Precompute per-term base powers
% And then take their product along the second axis
x_pow = power(x_row, p); % (m, n_x)
y_pow = power(y_row, q); % (m, n_x)
y_pow_prod = prod(y_pow, 2); % (m, 1)
x_pow_prod = prod(x_pow, 2); % (m, 1)
% Row-wise product excluding each column one by one.
x_pow_exc_prod = polynomial_bilevel.exclusive_row_products(x_pow);
y_pow_exc_prod = polynomial_bilevel.exclusive_row_products(y_pow);
switch with_respect_to
% Gradient with respect to x: (n_x, 1)
case 'x'
out = zeros(n_x, 1);
for k = 1:n_x
term = c .* y_pow_prod .* p(:,k) .* ...
x_pow_exc_prod(:,k) .* ( x_row(k) .^ max(p(:,k)-1, 0) );
out(k) = sum(term);
end
% Gradient with respect to y: (n_y, 1)
case 'y'
out = zeros(n_y, 1);
for j = 1:n_y
term = c .* x_pow_prod .* q(:,j) .* ...
y_pow_exc_prod(:,j) .* ( y_row(j) .^ max(q(:,j)-1, 0) );
out(j) = sum(term);
end
% Hessian with respect to x twice: (n_x, n_x)
case 'xx'
out = zeros(n_x, n_x);
for k = 1:n_x
for l = 1:n_x
if k == l
expk = p;
expk(:,k) = max(p(:,k) - 2, 0);
row_val = c .* y_pow_prod .* (p(:,k) .* max(p(:,k)-1, 0)) .* ...
prod(power(x_row, expk), 2);
else
expkl = p;
expkl(:,k) = max(p(:,k) - 1, 0);
expkl(:,l) = max(p(:,l) - 1, 0);
row_val = c .* y_pow_prod .* (p(:,k) .* p(:,l)) .* ...
prod(power(x_row, expkl), 2);
end
out(k,l) = sum(row_val);
end
end
% Hessian with respect to y twice. (n_x, n_y)
case 'yy'
out = zeros(n_y, n_y);
for j = 1:n_y
for h = 1:n_y
if j == h
expj = q;
expj(:,j) = max(q(:,j) - 2, 0);
row_val = c .* x_pow_prod .* (q(:,j) .* max(q(:,j)-1, 0)) .* ...
prod(power(y_row, expj), 2);
else
expjh = q;
expjh(:,j) = max(q(:,j) - 1, 0);
expjh(:,h) = max(q(:,h) - 1, 0);
row_val = c .* x_pow_prod .* (q(:,j) .* q(:,h)) .* ...
prod(power(y_row, expjh), 2);
end
out(j,h) = sum(row_val);
end
end
% Hessian with respect to x and y. (n_x, n_y)
case 'xy'
out = zeros(n_x, n_y);
for k = 1:n_x
X_part = p(:,k) .* x_pow_exc_prod(:,k) .* ( x_row(k) .^ max(p(:,k)-1, 0) );
for j = 1:n_y
Y_part = q(:,j) .* y_pow_exc_prod(:,j) .* ( y_row(j) .^ max(q(:,j)-1, 0) );
row_val = c .* (X_part .* Y_part);
out(k,j) = sum(row_val);
end
end
% Otherwise
otherwise
error('Unsupported derivative type: %s', with_respect_to);
end
end
function out = exclusive_row_products(matrix)
% Row-wise product excluding each column one by one.
%
% matrix: (m x n)
% out: (m x n), out(i,j) = prod_{k:k!=j} matrix(i,k)
%
[m, n] = size(matrix);
L = ones(m, n);
if n > 1
L(:, 2:end) = cumprod(matrix(:, 1:end-1), 2);
end
R = ones(m, n);
if n > 1
R(:, 1:end-1) = fliplr(cumprod(fliplr(matrix(:, 2:end)), 2));
end
out = L .* R;
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 JSON file
txt = fileread(path);
data = jsondecode(txt);
% Store the path. This is useful diagnostics
data.path = char(path);
% ================== REQUIRED keys ==================
% Objective coefficients are nessasary
% Convert the coefficents and exponents to doubles
keys = {'F', 'f'};
subkeys = {'c', 'p', 'q'};
for i = 1:length(keys)
key = keys{i};
% Diagnostics message
msg = sprintf(strcat('Invalid polynomial for %s. ', ...
'It should look like this ', ...
'\n "%s": { "c": [[]], "p": [[]], "q": [[]] }. ', ...
'\n See JSON file: %s.'), key, key, path ...
);
assert(isfield(data, key), msg)
for j = 1:length(subkeys)
subkey = subkeys{j};
assert(isfield(data.(key), subkey), msg)
data.(key).(subkey) = double(data.(key).(subkey));
end
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=- %
%% 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_map = dictionary( ...
'x', data.n_x, ...
'y', data.n_y, ...
'F', 1, ...
'G', length(data.ineq_upper) + length(data.rhs_upper) + length(data.x_upper_bound) + length(data.x_lower_bound), ...
'H', 0, ...
'f', 1, ...
'g', length(data.ineq_lower) + length(data.rhs_lower) + length(data.y_upper_bound) + length(data.y_lower_bound), ...
'h', 0 ...
);
if isKey(n_map, key)
n = n_map(key);
else
n = n_map;
end
end
end
endimport json
import os.path
from bolib3 import np
"""
Polynomial Bilevel Program (PBP)
================================
The objective function and constraints are polynomials.
Each polynomial is parameterised by m, c, p, q as defined below.
These parameters can be loaded from the JSON files in bolib3/data/polynomial_bilevel.
Polynomial expression
=====================
m: Number of terms.
c: Vector (m,) of coefficients.
p: Matrix (m, n_x) of exponents for x.
q: Matrix (m, n_y) of exponents for y.
In LaTeX syntax:
Polynomial(x, y) := sum_{t=1,...,m}(
c_t
* product_{i=1,...,n_x}( x_i^{p_{t,i}} )
* product_{i=1,...,n_y}( y_i^{q_{t,i}} )
)
In python syntax:
polynomial(x, y) := term_1 + term_2 + ... + term_m.
where term_i := c[i] * (x[0]**p[i,0] * x[1]**p[i,1] * ...) * (y[0]**q[i,0] * y[1]**q[i,1] * ...)
"""
# Properties
name: str = "polynomial_bilevel"
category: str = "archetype"
subcategory: str = ""
author: str = "Samuel Ward"
# Methods
def F(x, y, data=None):
"""Upper-level objective function"""
return polynomial(x, y, **data['F'])
def G(x, y, data=None):
"""Upper-level inequality constraints"""
ineq_list = [np.empty(0)]
# Polynomial constraints
if len(data["ineq_upper"]) > 0:
ineq_list.append(np.array([
polynomial(x, y, **q) for q in data['ineq_upper']
]))
# Linear constraints
if data["rhs_upper"].size > 0:
ineq_list.append((data["A_upper"]@x) + (data["B_upper"]@y) - data["rhs_upper"])
# Lower bounds on x
if data["x_lower_bound"].size > 0:
ineq_list.append(x - data["x_lower_bound"])
# Upper bounds on x
if data["x_upper_bound"].size > 0:
ineq_list.append(data["x_upper_bound"] - x)
# All 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=None):
"""Lower-level objective function"""
return polynomial(x, y, **data['f'])
def g(x, y, data=None):
"""lower-level inequality constraints"""
ineq_list = [np.empty(0)]
# Polynomial constraints
if len(data["ineq_lower"]) > 0:
ineq_list.append(np.array([
polynomial(x, y, **q) for q in data['ineq_lower']
]))
# Linear constraints
if data["rhs_lower"].size > 0:
ineq_list.append((data["A_lower"]@x) + (data["B_lower"]@y) - data["rhs_lower"])
# Lower bounds on y
if data["y_lower_bound"].size > 0:
ineq_list.append(y - data["y_lower_bound"])
# Upper bounds on y
if data["y_upper_bound"].size > 0:
ineq_list.append(data["y_upper_bound"] - y)
# All constraints G(x, y) >= 0
return np.concatenate(ineq_list)
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=-
# Polynomial Expression
# -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-
# =============================================================
def polynomial(x, y, c, p, q, **kwargs):
"""
Evaluates the polynomial expression
sum_{t=1,...,m} ( c_t
* product_{i=1,...,n_x}( x_i^{p_{t,i}} )
* product_{i=1,...,n_y}( y_i^{q_{t,i}} )
)
:param x: Array of values for decision variable x of shape (n_x)
:param y: Array of values for decision variable y of shape (n_y)
:param c: Coefficients of shape (m, )
:param p: Exponents of x shape (m, n_x)
:param q: Exponents of y shape (m, n_y)
"""
return np.sum(c*np.prod(np.power(x, p), axis=-1)*np.prod(np.power(y, q), axis=-1))
def _d_polynomial_dx(x, y, c, p, q, **kwargs):
"""First-order derivative of the polynomial (see above) with respect to x"""
m, n_x = p.shape
gradient_x = np.zeros(n_x)
y_term = np.prod(np.power(y, q), axis=-1)
for i in range(n_x):
not_i = ~np.isin(range(n_x), i)
x_term = p[:, i]*np.power(x[i], p[:, i] - 1)*np.prod(np.power(x[not_i], p[:, not_i]), axis=-1)
gradient_x[i] = np.sum(c*x_term*y_term)
return gradient_x
def _d_polynomial_dy(x, y, c, p, q, **kwargs):
"""First-order derivative of the polynomial (see above) with respect to y"""
m, n_y = q.shape
gradient_y = np.zeros(n_y)
x_term = np.prod(np.power(x, p), axis=-1)
for i in range(n_y):
not_i = ~np.isin(range(n_y), i)
y_term = q[:, i]*np.power(y[i], q[:, i] - 1)*np.prod(np.power(y[not_i], q[:, not_i]), axis=-1)
gradient_y[i] = np.sum(c*x_term*y_term)
return gradient_y
def _d2_polynomial_dxx(x, y, c, p, q, **kwargs):
"""Second-order derivative (Hessian) of the polynomial (see above) with respect to x twice"""
m, n_x = p.shape
hessian_xx = np.zeros((n_x, n_x))
y_term = np.prod(np.power(y, q), axis=1)
for i in range(n_x):
for j in range(n_x):
not_ij = ~np.isin(range(n_x), (i, j))
x_term = np.prod(np.power(x[not_ij], p[:, not_ij]), axis=-1)
if i == j:
x_term *= p[:, i]*(p[:, i] - 1)*np.power(x[i], p[:, i] - 2)
else:
x_term *= p[:, i]*p[:, j]*np.power(x[i], p[:, i] - 1)*np.power(x[j], p[:, j] - 1)
hessian_xx[i, j] = np.sum(c*y_term*x_term)
return hessian_xx
def _d2_polynomial_dyy(x, y, c, p, q, **kwargs):
"""Second-order derivative (Hessian) of the polynomial (see above) with respect to y twice"""
m, n_y = q.shape
hessian_yy = np.zeros((n_y, n_y))
x_term = np.prod(np.power(x, p), axis=1)
for i in range(n_y):
for j in range(n_y):
not_ij = ~np.isin(range(n_y), (i, j))
y_term = np.prod(np.power(y[not_ij], q[:, not_ij]), axis=-1)
if i == j:
y_term *= q[:, i]*(q[:, i] - 1)*np.power(y[i], q[:, i] - 2)
else:
y_term *= q[:, i]*np.power(y[i], q[:, i] - 1)*q[:, j]*np.power(y[j], q[:, j] - 1)
hessian_yy[i, j] = np.sum(c*y_term*x_term)
return hessian_yy
def _d2_polynomial_dxy(x, y, c, p, q, **kwargs):
"""Second-order derivative (Hessian) of the polynomial (see above) with respect to x and y"""
m, n_x = p.shape
m, n_y = q.shape
hessian_xy = np.zeros((n_x, n_y))
for i in range(n_x):
for j in range(n_y):
term = c*p[:, i]*np.power(x[i], p[:, i] - 1)*q[:, j]*np.power(y[j], q[:, j] - 1)
for k in range(n_x):
if k != i:
term *= np.power(x[k], p[:, k])
for l in range(n_y):
if l != j:
term *= np.power(y[l], q[:, l])
hessian_xy[i, j] = np.sum(term)
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# -=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 _d_polynomial_dx(x, y, **data['F'])
def dF_dy(x, y, data):
"""First-order derivative of the upper-level objective function F(x,y) with respect to y"""
return _d_polynomial_dy(x, y, **data['F'])
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 _d2_polynomial_dxx(x, y, **data['F'])
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 _d2_polynomial_dyy(x, y, **data['F'])
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 _d2_polynomial_dxy(x, y, **data['F'])
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 len(data["ineq_upper"]):
jacobian_list.append(np.array([
_d_polynomial_dx(x, y, **p) for p in data['ineq_upper']
]))
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 len(data["ineq_upper"]):
jacobian_list.append(np.array([
_d_polynomial_dy(x, y, **p) for p in data['ineq_upper']
]))
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 _d_polynomial_dx(x, y, **data['f'])
def df_dy(x, y, data):
"""First-order derivative of the lower-level objective function f(x,y) with respect to y"""
return _d_polynomial_dy(x, y, **data['f'])
def d2f_dxx(x, y, data):
"""Second-order derivative (Hessian) of the lower-level objective function f(x,y) with respect x twice"""
return _d2_polynomial_dxx(x, y, **data['f'])
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 _d2_polynomial_dyy(x, y, **data['f'])
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 _d2_polynomial_dxy(x, y, **data['f'])
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 len(data["ineq_lower"]):
jacobian_list.append(np.array([
_d_polynomial_dx(x, y, **p) for p in data['ineq_lower']
]))
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 len(data["ineq_lower"]):
jacobian_list.append(np.array([
_d_polynomial_dy(x, y, **p) for p in data['ineq_lower']
]))
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=-
# =============================================================
datasets: list = sorted([
'dempe_dutta_2012_ex24.json', 'colson_2002_bipa1.json', 'colson_2002_bipa2.json',
'colson_2002_bipa4.json', 'floudas_zlobec_1998.json', 'gumus_floudas_2001_ex1.json',
'ishizuka_aiyoshi_1992a.json', 'kleniati_adjiman_2014_ex3.json', 'lampariello_sagratella_2017_ex23.json',
'mitsos_barton_2006_ex39.json', 'mitsos_barton_2006_ex310.json', 'mitsos_barton_2006_ex311.json',
'mitsos_barton_2006_ex312.json', 'mitsos_barton_2006_ex313.json', 'mitsos_barton_2006_ex314.json',
'mitsos_barton_2006_ex315.json', 'mitsos_barton_2006_ex316.json', 'mitsos_barton_2006_ex317.json',
'mitsos_barton_2006_ex318.json', 'mitsos_barton_2006_ex319.json', 'mitsos_barton_2006_ex320.json',
'mitsos_barton_2006_ex321.json', 'mitsos_barton_2006_ex322.json', 'mitsos_barton_2006_ex323.json',
'mitsos_barton_2006_ex324.json', 'mitsos_barton_2006_ex325.json', 'mitsos_barton_2006_ex326.json',
'nie_wang_ye_2017_ex34.json', 'nie_wang_ye_2017_ex52.json', 'nie_wang_ye_2017_ex54.json',
'nie_wang_ye_2017_ex57.json', 'nie_wang_ye_2017_ex58.json', 'nie_wang_ye_2017_ex61.json',
'outrata_1993_ex31.json', 'outrata_1993_ex32.json', 'outrata_1994_ex31.json',
'paulavicius_adjiman_2017a.json', 'paulavicius_adjiman_2017b.json', 'shimizu_etal_1997b.json',
'vogel_2012.json', 'ye_zhu_2010_ex42.json', 'ye_zhu_2010_ex43.json',
])
paths: list = [
os.path.join("bolib3", "data", "polynomial_bilevel", dataset) for dataset in datasets
]
def read_data(path: str = datasets[0], verbose: int = 0):
"""
Must give a path to a JSON file that follows the format of 'polynomial_template.json'.
For example, it should begin something like this:
{
"name": "author_1990",
"n_x": 3,
"n_y": 3,
"F": {
"m": 2,
"c": [1.0, 2.0],
"p": [[11, 12, 13], [21, 22, 23]],
"q": [[11, 12, 13], [21, 22, 23]],
},
"""
# 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
for key in ('n_x', 'n_y', 'F', 'f'):
assert key in data, f"Missing required key \"{key}\". See JSON file: {path}."
# Cast all the polynomial parameters to numpy arrays
for key1, poly in list_polynomials(data):
for key2 in ('c', 'p', 'q'):
assert key2 in poly, (f"Missing required subkey \"{key1}\"{{ \"{key2}\": <value>, ... }}. "
f"See JSON file: {path}.")
poly[key2] = np.array(poly[key2])
# ================== 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)
# Validate the dimensions of all arrays
_validate_data_dimensions(data)
# Print out the program
if verbose:
display_program(data)
return data
def _validate_data_dimensions(data: dict):
"""
Checks the shapes of all the matrices and vectors in the data are consistent.
"""
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', '')))
# Check that the dimensions of the matrices for the linear constraints are consistent
for key, expected_shape in {
'A_upper': (n_G, n_x),
'B_upper': (n_G, n_y),
'rhs_upper': (n_G,),
'A_lower': (n_g, n_x),
'B_lower': (n_g, 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}'."
# Check that the terms of the polynomials are consistent
for key1, poly in list_polynomials(data):
m = poly['m']
for key2, expected_shape in {
'c': (m,),
'p': (m, n_x),
'q': (m, n_y)
}.items():
observed_shape = poly[key2].shape
assert observed_shape == expected_shape, \
(f"Dimension of {key1}[\"{key2}\"] is {observed_shape} "
f"but was expect to be {expected_shape}. See: '{identifier}'.")
def list_polynomials(data: dict):
"""Returns a list of (key:str, polynomial:dict) tuples"""
for key in ('F', 'f', 'ineq_upper', 'ineq_lower'):
assert key in data, f"Missing required key \"{key}\"."
return [
('F', data['F'])] + [
('f', data['f'])] + [
(f'ineq_upper[{i}]', data['ineq_upper'][i]) for i in range(len(data['ineq_upper']))] + [
(f'ineq_lower[{i}]', data['ineq_lower'][i]) for i in range(len(data['ineq_lower']))
]
def display_polynomial(c, p, q, index_from_1=1, **kwargs):
"""
Converts a polynomial to a markdown/LaTeX textual representation.
e.g. "16.0 x_1 y_1^4 + 2.0 x_1 y_1^3 -8.0 x_1 y_1^2 -1.5 x_1 y_1 +0.5 x_1"
:param c: Coefficients of shape (m, )
:param p: Exponents of x shape (m, n_x)
:param q: Exponents of y shape (m, n_y)
:param index_from_1: Set to 1 for one-based indexing like MATLAB or to 0 for zero-based indexing like python.
"""
text = ""
# Iterate terms
for term in range(len(c)):
text += f"{float(c[term]):>+8.2f}"
# Iterate variables
variables = []
for v, exponents_array in {'x': p[term], 'y': q[term]}.items():
for i, exponent in enumerate(exponents_array):
index = f"_{i + index_from_1}" if len(exponents_array) > 1 else ""
if exponent == 1:
variables.append(f"{v}{index}")
elif exponent != 0:
variables.append(f"{v}{index}^{exponent}")
if len(variables):
text += ' (' + ' '.join(variables) + ') '
return text
def display_constraints(data: dict, key: str = "ineq_upper"):
constraints = [display_polynomial(**data[key][i]) for i in range(len(data[key]))]
width = max([0.0] + [len(c) for c in constraints])
return ''.join(f"\n{' '*12}{c:<{width}} >= 0.000" for c in constraints)
def display_program(data):
return f"""
Upper-level Program
===================
minimise {display_polynomial(**data["F"])}
subject to Ax + By >= rhs{display_constraints(data, "ineq_upper")}
{data["x_lower_bound"]} <= x <= {data["x_upper_bound"]}
Lower-level Program
===================
minimise {display_polynomial(**data["f"])}
subject to Ax + By >= rhs{display_constraints(data, "ineq_lower")}
{data["y_lower_bound"]} <= y <= {data["y_upper_bound"]}
"""
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# -=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['ineq_upper']) + len(data['rhs_upper']) + len(data['x_upper_bound']) + len(data['x_lower_bound']),
"H": 0,
"f": 1,
"g": len(data['ineq_lower']) + 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