Added 30/09/2025
Archetype
quadratically_constrained_quadratic_bilevel
Datasets
Description
Allende, Gemayqzel Bouza and Still, Georg (2013) (see page 331).
Solving bilevel programs with the KKT-approach.
https://doi.org/10.1007/s10107-012-0535-xDimension
{
"x": 2,
"y": 2,
"F": 1,
"G": 5,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "global",
"x": [0.5,0.5],
"y": [0.5,0.5],
"F": -1.5,
"G": [1.5,0.5,0.5,0.5,0.5],
"H": [],
"f": -0.5,
"g": [0,0],
"h": []
}Description
Bard, Jonathan F. (1988) (see page 24, Example 3).
Convex two-level optimization.
https://doi.org/10.1007/BF01580720Dimension
{
"x": 2,
"y": 2,
"F": 1,
"G": 3,
"H": 0,
"f": 1,
"g": 4,
"h": 0
}Solution
{
"optimality": "unknown",
"x": [1.45,0.95],
"y": [1.88,0.64],
"F": -12.68,
"G": [-0.0024999999999995026,1.45,0.95],
"H": [],
"f": -1.02,
"g": [-0.01499999999999968,0.02999999999999936,1.88,0.64],
"h": []
}Description
Dempe, S. (1992).
A necessary and a sufficient optimality condition for bilevel programming problems.
https://doi.org/10.1080/02331939208843831Dimension
{
"x": 2,
"y": 2,
"F": 1,
"G": 1,
"H": 0,
"f": 1,
"g": 2,
"h": 0
}Solution
{
"optimality": "best_known",
"x": [0,0],
"y": [0,-0.5],
"F": -0.5,
"G": [0],
"H": [],
"f": 0.625,
"g": [0,0],
"h": []
}Description
Dempe, S. (1992).
A necessary and a sufficient optimality condition for bilevel programming problems.
https://doi.org/10.1080/02331939208843831Dimension
{
"x": 1,
"y": 1,
"F": 1,
"G": 0,
"H": 0,
"f": 1,
"g": 1,
"h": 0
}Solution
{
"optimality": "unknown",
"x": [0],
"y": [0],
"F": 31.25,
"G": [],
"H": [],
"f": 4,
"g": [0],
"h": []
}\subsection{quadratically\_constrained\_quadratic\_bilevel}
\label{subsec:quadratically_constrained_quadratic_bilevel}
% Quadratic Expression
Let $\x\in\R^{n_x},\y\in\R^{n_y}$ be the lower and upper level decision variables, respectively.
Let
$c\in\R^{n_x}$,
$d\in\R^{n_y}$
be data vectors and
$P\in\R^{n_x\times n_x}$,
$Q\in\R^{n_y\times n_y}$,
$R\in\R^{n_x\times n_y}$
data matrices and
$s\in\R$ a scalar
then the \emph{Quadratic Expression} in $x,y$ is defined as:
\begin{equation*}
\mathcal{Q}_{c,d,P,Q,R,s}(x,y)\coloneqq
\frac{1}{2}
\begin{bmatrix}x\\y\end{bmatrix}^\top
\begin{bmatrix}
P & R\\
R^\top & Q
\end{bmatrix}
\begin{bmatrix}x\\y\end{bmatrix}
+
\begin{bmatrix}c\\d\end{bmatrix}^\top
\begin{bmatrix}x\\y\end{bmatrix}
+s.
\end{equation*}
% Quadratically Constrained Quadratic Bilevel Program
Given a sequence of quadratic expression $\mathcal{Q}^0_\upper, \mathcal{Q}^k_\upper, \mathcal{Q}^0_\low, \mathcal{Q}^k_\low$ the \emph{Quadratically Constrained Quadratic Bilevel Program}~\eqref{eq:QCQBP} over $x,y$ is defined as:
\begin{flalign*}
\label{eq:QCQBP}
\tag{QCQBP}
\minimise_{\x\in\R^{n_x}, \y\in\R^{n_y}} \quad
& \mathcal{Q}^0_\upper(x,y)\\
\subjectto \quad
& \mathcal{Q}^k_\upper(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}
\mathcal{Q}^0_{\low}(x,y)\\
\st[1]\quad\mathcal{Q}^k_{\low}(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*}
% Example
We store the parameters of each quadratic expression ($c, d, P, Q, R, s$) as named vectors and matrices in a JSON file.
For example consider the~\eqref{eq:QCQBP} instance given by Allende and Still~\cite[page~331]{Allende2013}:
\begin{flalign*}
\minimise_{\x\in\R^{n_x}, \y\in\R^{n_y}} \quad
& F(x,y)\coloneqq -2x_1 -2x_2 -x_1^2 +x_2^2 +y_1^2+y_2^2\\
\subjectto \quad
& x\geq0,\quad y\geq 0,\quad -x_1 + 2 \geq 0,\\
& y \in \argmin_{y\in\R^{n_y}}
\left\{
\begin{array}{l}
y_1^2 - 2x_1y_1 + y_2^2 - 2x_2y_2\\
\st[1]\quad 0.25 - (y_1 - 1)^2 \;\geq 0,\\
\st[2]\quad 0.25 - (y_2 - 1)^2 \;\geq 0.\\
\end{array}
\right.
\end{flalign*}
% JSON
The upper-level objective is a quadratic expression whose coefficients are stored in the following way:
\begin{lstlisting}[style=custom,language=JSON,label={lst:lstlisting}]
"F": {
"c": [-2, -2],
"d": [0, 0],
"P": [[-2, 0], [0, 2]],
"Q": [[2, 0], [0, 2]],
"R": [[0, 0], [0, 0]],
"s": 0
},
\end{lstlisting}classdef quadratically_constrained_quadratic_bilevel
% Quadratically Constrained Quadratic Bilevel Program (QCQBP)
% ===========================================================
% The objective function and constraints are quadratic expression in (x, y).
% Each quadratic expression is parameterised by (c, d, P, R, Q, s).
%
%
% Quadratic expression in x and y
% ===============================
% x: Leader's decision variables; vector of shape (n_x,).
% y: Follower's decision variables; vector of shape (n_y,).
% c: Linear coefficient of x. vector of shape (n_x,).
% d: Linear coefficient of y. shape (n_y,).
% P: Quadratic form matrix of shape (n_x, n_x).
% Q: Quadratic form matrix of shape (n_y, n_y).
% R: Quadratic form matrix of shape (n_x, n_y).
% s: Scalar constant.
%
%
% Q(x,y) := 0.5 [x]'* [ P R ] * [x] + [c]'[x] + s
% [y] [ R' Q ] [y] [d] [y]
%
properties(Constant)
name = 'quadratically_constrained_quadratic_bilevel';
category = 'archetype';
subcategory = '';
author = 'Samuel Ward';
datasets = {'allende_still_2013.json', 'bard_1988_ex3.json', 'dempe_1992a.json', 'dempe_1992b.json'};
paths = fullfile('bolib3', 'data', 'quadratically_constrained_quadratic_bilevel', quadratic_bilevel.datasets);
end
methods(Static)
function val = quadratic(x, y, coefficients)
% Quadratic in
% x: leader's decision variables; vector of shape (n_x,).
% y: follower's decision variables; vector of shape (n_y,).
%
% The coefficients argument should be a stuct with feilds
% c: linear coefficient of x. vector of shape (n_x,).
% d: linear coefficient of y. shape (n_y,).
% P: quadratic form matrix of shape (n_x, n_x).
% Q: quadratic form matrix of shape (n_y, n_y).
% R: quadratic form matrix of shape (n_x, n_y).
% s: scalar constant.
val = 0.5 * x' * coefficients.P * x + ...
0.5 * y' * coefficients.Q * y + ...
1.0 * x' * coefficients.R * y + ...
dot(coefficients.c, x) + ...
dot(coefficients.d, y) + ...
coefficients.s;
end
function val = F(x, y, data)
% Upper-level objective function
val = quadratically_constrained_quadratic_bilevel.quadratic( ...
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) = quadratically_constrained_quadratic_bilevel.quadratic( ...
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 = quadratically_constrained_quadratic_bilevel.quadratic( ...
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) = quadratically_constrained_quadratic_bilevel.quadratic( ...
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
% Dimention (n_x × 1)
dx = (data.F.P * x) + (data.F.R * y) + data.F.c;
end
function dy = dF_dy(x, y, data)
% ∂F/∂y
% Dimention (n_y × 1)
dy = (data.F.Q * y) + (data.F.R' * x) + data.F.d;
end
function Hxx = d2F_dxx(~, ~, data)
% ∂²F/∂x²
% Dimention (n_x × n_x)
Hxx = data.F.P;
end
function Hyy = d2F_dyy(~, ~, data)
% ∂²F/∂y²
% Dimention (n_y × n_y)
Hyy = data.F.Q;
end
function Hxy = d2F_dxy(~, ~, data)
% ∂²F/∂x∂y
% Dimention (n_x × n_y)
Hxy = data.F.R;
end
function Jx = dG_dx(x, y, data)
% Jacobian of G wrt x
% Dimention (|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)
ineq = data.ineq_upper(i);
Jx(i,:) = (ineq.P * x) + (ineq.R * y) + ineq.c;
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:
% Dimention (|G| × n_y)
Jy = zeros(0, data.n_y);
if ~isempty(data.ineq_upper)
Jy = zeros(length(data.ineq_upper), data.n_x);
for i = 1:length(data.ineq_upper)
ineq = data.ineq_upper(i);
Jy(i,:) = (ineq.Q * y) + (ineq.R' * x) + ineq.d;
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
% 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(x, y, data)
% ∂f/∂x
% Dimention (n_x × 1)
dx = (data.f.P * x) + (data.f.R * y) + data.f.c;
end
function dy = df_dy(x, y, data)
% ∂f/∂y
% Dimention (n_y × 1)
dy = (data.f.Q * y) + (data.f.R' * x) + data.f.d;
end
function Hxx = d2f_dxx(~, ~, data)
% ∂²f/∂x²
% Dimention (n_x × n_x)
Hxx = data.f.P;
end
function Hyy = d2f_dyy(~, ~, data)
% ∂²f/∂y²
% Dimention (n_y × n_y)
Hyy = data.f.Q;
end
function Hxy = d2f_dxy(~, ~, data)
% ∂²f/∂x∂y
% Dimention (n_x × n_y)
Hxy = data.f.R;
end
function Jx = dg_dx(x, y, data)
% Jacobian of g wrt x
% Dimention (|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)
ineq = data.ineq_lower(i);
Jx(i,:) = (ineq.P * x) + (ineq.R * y) + ineq.c;
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
% Dimention (|g| × n_y)
Jy = zeros(0, data.n_y);
if ~isempty(data.ineq_lower)
Jy = zeros(length(data.ineq_lower), data.n_x);
for i = 1:length(data.ineq_lower)
ineq = data.ineq_lower(i);
Jy(i,:) = (ineq.Q * y) + (ineq.R' * x) + ineq.d;
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
% 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 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 to double data type
keys = {'F', 'f'};
subkeys = {'c', 'd', 'P', 'Q', 'R', 's'};
for i = 1:length(keys)
key = keys{i};
% Diagnostics message
msg = sprintf(strcat('Invalid quadratic form for %s. ', ...
'It should look like this ', ...
'\n "%s": { "c": [], "d": [], "P": [[]], "Q": [[]], "R": [[]], "s": 0 } ', ...
'\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 = 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,key)
n = n(key);
end
end
end
endimport json
import os
from bolib3 import np
"""
Quadratically Constrained Quadratic Bilevel Program (QCQBP)
===========================================================
The objective function and constraints are quadratic expression in (x, y).
Each quadratic expression is parameterised by (c, d, P, R, Q, s).
Quadratic expression in x and y
===============================
x: Leader's decision variables; vector of shape (n_x,).
y: Follower's decision variables; vector of shape (n_y,).
c: Linear coefficient of x. vector of shape (n_x,).
d: Linear coefficient of y. shape (n_y,).
P: Quadratic form matrix of shape (n_x, n_x).
Q: Quadratic form matrix of shape (n_y, n_y).
R: Quadratic form matrix of shape (n_x, n_y).
s: Scalar constant.
Q(x,y) := 0.5 [x]' [ P R ] [x] + [c]'[x] + s
[y] [ R' Q ] [y] [d] [y]
Q(x,y) := 0.5*(x'Px) + 0.5*(y'Qy) + (x'Ry) + c'x + d'y + s
"""
# Properties
name: str = "quadratically_constrained_quadratic_bilevel"
category: str = "archetype"
subcategory: str = ""
author: str = "Samuel Ward"
# Methods
def F(x, y, data=None):
"""Upper-level objective function"""
return quadratic(x, y, **data['F'])
def G(x, y, data=None):
"""Upper-level inequality constraints"""
ineq_list = [np.empty(0)]
# Quadratic constraints
if len(data["ineq_upper"]):
ineq_list.append(np.array([
quadratic(x, y, **q) for q in data['ineq_upper']
]))
# Linear constraints
if data["rhs_upper"].size:
ineq_list.append(data["A_upper"]@x + data["B_upper"]@y - data["rhs_upper"])
# Lower bounds on x
if data["x_lower_bound"].size:
ineq_list.append(x - data["x_lower_bound"])
# Upper bounds on x
if data["x_upper_bound"].size:
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 quadratic(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"]):
ineq_list.append(np.array([
quadratic(x, y, **q) for q in data['ineq_lower']
]))
# Linear constraints
if data["rhs_lower"].size:
ineq_list.append(data["A_lower"]@x + data["B_lower"]@y - data["rhs_lower"])
# Lower bounds on y
if data["y_lower_bound"].size:
ineq_list.append(y - data["y_lower_bound"])
# Upper bounds on y
if data["y_upper_bound"].size:
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=-
# Quadratic Expression
# -=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-=x=-
# =============================================================
def quadratic(x, y, c, d, P, Q, R, s, **kwargs):
"""
Evaluates the quadratic expression.
Q(x, y) := 0.5*(x' P x) + 0.5*(y' Q y) + 1.0*(x' R y) + c'x + d'y + s
: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: Linear coefficient of x. shape (n_x,).
:param d: Linear coefficient of y. shape (n_y,).
:param P: Quadratic form matrix of shape (n_x, n_x).
:param Q: Quadratic form matrix of shape (n_y, n_y).
:param R: Quadratic form matrix of shape (n_x, n_y).
:param s: Scalar constant.
"""
return np.sum((
0.5*x@P@x,
0.5*y@Q@y,
1.0*x@R@y,
np.dot(c, x),
np.dot(d, y),
s
))
def _d_quadratic_dx(x, y, c, d, P, Q, R, s, **kwargs):
"""Derivative of quadratic (see above) with respect to x"""
return P@x + R@y + c
def _d_quadratic_dy(x, y, c, d, P, Q, R, s, **kwargs):
"""Derivative of quadratic (see above) with respect to y"""
return Q@y + x@R + d
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# -=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_quadratic_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_quadratic_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 data['F']['P']
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 data['F']['Q']
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 data['F']['R']
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_quadratic_dx(x, y, **q) for q 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_quadratic_dy(x, y, **q) for q 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 data['f']['P']@x + data['f']['R']@y + data['f']['c']
def df_dy(x, y, data):
"""First-order derivative of the lower-level objective function f(x,y) with respect to y"""
return data['f']['Q']@y + x@data['f']['R'] + data['f']['d']
def d2f_dxx(x, y, data):
"""Second-order derivative (Hessian) of the lower-level objective function f(x,y) with respect x twice"""
return data['f']['P']
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 data['f']['Q']
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 data['f']['R']
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_quadratic_dx(x, y, **q) for q 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_quadratic_dy(x, y, **q) for q 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 = [
'allende_still_2013.json', 'bard_1988_ex3.json', 'dempe_1992a.json', 'dempe_1992b.json'
]
paths: list = [
os.path.join("bolib3", "data", "quadratically_constrained_quadratic_bilevel", d) for d in datasets
]
def read_data(path=''):
"""
Must give a path to a JSON file that follows the format of quadratically_constrained_quadratic_bilevel.json.
Keys "F" and "f" should have sub-keys "c", "d", "P", "Q", "R", "s" defining the coefficients of a quadratic.
For example, it should begin something like this:
{
"name": "author_1990",
"n_x": 2,
"n_y": 2,
"F": {
"c": [1.0, 2.0],
"d": [1.0, 2.0],
"P": [[11.0, 12.0], [21.0, 22.0]],
"Q": [[11.0, 12.0], [21.0, 22.0]],
"R": [[11.0, 12.0], [21.0, 22.0]],
"s": 0.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
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 quadratic expression parameters to numpy arrays
for key1, quad in list_quadratic_expressions(data):
for key2 in ("c", "d", "P", "Q", "R", "s"):
assert key2 in quad, (f"Missing required subkey \"{key1}\"{{ \"{key2}\": <value>, ... }}. "
f"See JSON file: {path}.")
quad[key2] = np.array(quad[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)
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_quadratic_expressions(data):
for key2, expected_shape in {
"c": (n_x,),
"d": (n_y,),
"P": (n_x, n_x),
"Q": (n_y, n_y),
"R": (n_x, n_y),
"s": ()
}.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_quadratic_expressions(data: dict):
"""Returns a list of (key:str, quadratic: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']))
]
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# -=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