Bilevel via KKT + Model.complementarity()#
A bilevel (Stackelberg) program has a leader optimizing over a follower that solves its own problem in response:
When the lower level is convex in \(y\) (for each fixed \(x\)), the follower’s
KKT conditions are necessary and sufficient for its optimality, so
y ∈ argmin{…} can be replaced by a single-level system
[Bard, 1998, Colson et al., 2007, Dempe, 2002]:
stationarity \(\nabla_y L = 0\), with \(L = f + \sum_i \mu_i\, g_i\),
primal feasibility \(g_i(x,y) \le 0\),
dual feasibility \(\mu_i \ge 0\),
complementary slackness \(0 \le \mu_i \perp -g_i(x,y) \ge 0\).
The result is a Mathematical Program with Equilibrium Constraints (MPEC) [Luo et al., 1996]. Its hard part is the last line: the complementary-slackness product \(\mu_i\,g_i = 0\) is nonconvex, and the classic remedy is the Fortuny-Amat–McCarl big-M disjunction [Fortuny-Amat and McCarl, 1981].
This notebook wires that last step through discopt’s fluent
Model.complementarity(), which encodes each pair as the exact disjunction
\([\mu_i = 0] \vee [g_i = 0]\) and lets the GDP machinery lower it
[Raman and Grossmann, 1994]. The companion notebooks cover the pieces
separately: complementarity & MPEC for the
complementarity() API on standalone MPCCs, and bilevel
optimization for the automated BilevelProblem reduction. Here
we tie them together end-to-end: bilevel → KKT → complementarity() on the
complementary-slackness pairs, to a certified global optimum.
import numpy as np
import discopt.modeling as dm
from discopt.bilevel import diff # symbolic derivative over the expression DAG
def value(res, name):
"""Scalar solution value for a variable by name."""
return float(np.asarray(res.x[name]))
A small Stackelberg problem#
A regulator (leader) sets an allotment \(x \in [0, 5]\) and wants the pair \((x, y)\) close to a target — objective \(F = (x-3)^2 + (y-1)^2\). A firm (follower) then picks its activity level \(y\) to minimise a convex cost, but may not exceed the allotment (\(y \le x\)):
The follower’s unconstrained optimum is \(y = 4\), so the coupling cap \(y \le x\) binds whenever \(x < 4\). Substituting the follower response \(y = \min(4, x)\) into \(F\), the leader’s best choice is \(x^\star = 2\), giving \(y^\star = 2\) and \(F^\star = 2\) — with the cap \(y \le x\) active at the optimum (a multiplier \(\mu^\star = 2 > 0\)). Because the cap is active, complementary slackness is not vacuous here: it is exactly what pins the solution.
def build_kkt_model():
# Leader + follower KKT reduction, up to (not including) complementary slackness.
# Emits stationarity and primal feasibility onto the model and returns the
# complementary-slackness ingredients (mu, g) for the caller to encode.
m = dm.Model("stackelberg")
x = m.continuous("x", lb=0, ub=5) # leader: allotment
y = m.continuous("y", lb=0, ub=5) # follower: activity level
m.minimize((x - 3) ** 2 + (y - 1) ** 2) # leader objective F(x, y)
# Follower (convex QP in y): min 0.5 y^2 - 4 y s.t. g = y - x <= 0, y >= 0.
f_lower = 0.5 * y * y - 4 * y
g = y - x # the coupling constraint written as g <= 0
# Dual multiplier for g. A KKT multiplier has no a-priori upper bound; a finite
# bound is a modeling assertion (like any big-M) that must be >= the true dual.
# Here mu* = 2, so ub=20 is safely valid and keeps the disjunction's big-M finite.
mu = m.continuous("mu", lb=0.0, ub=20.0)
# Stationarity: dL/dy = 0 with L = f_lower + mu * g. `diff` returns an ordinary
# expression, so this stays on discopt's certified global path.
lagrangian = f_lower + mu * g
m.subject_to(diff(lagrangian, y) == 0, name="stationarity")
# Primal feasibility of the follower constraint.
m.subject_to(g <= 0, name="primal_feas")
return m, x, y, mu, g
# Inspect the emitted KKT conditions (pure modeling — no solve yet).
_m, _x, _y, _mu, _g = build_kkt_model()
print("stationarity dL/dy =", diff(0.5 * _y * _y - 4 * _y + _mu * _g, _y), "== 0")
print("primal feas. g = y - x <= 0")
print("dual feas. mu >= 0")
print("compl. slackness 0 <= mu _|_ (x - y) >= 0")
stationarity dL/dy = ((((0.5 * y) + (0.5 * y)) - 4) + mu) == 0
primal feas. g = y - x <= 0
dual feas. mu >= 0
compl. slackness 0 <= mu _|_ (x - y) >= 0
Encode complementary slackness with Model.complementarity()#
Everything above is smooth and convex. The only equilibrium condition left is
the complementary slackness \(0 \le \mu \perp -g \ge 0\), i.e. \(0 \le \mu \perp
(x - y) \ge 0\). We hand it to Model.complementarity(mu, -g), whose default
method="gdp" encodes it as the exact disjunction \([\mu = 0] \vee [-g = 0]\) and
lowers it to a big-M MILP with a selector binary. Branch-and-bound then branches
on the finite either/or — the follower is either off the cap (\(\mu = 0\))
or on it (\(g = 0\)) — and certifies the result.
m, x, y, mu, g = build_kkt_model()
m.complementarity(mu, -g, name="compl_slack") # 0 <= mu _|_ (x - y) >= 0, default gdp
res_gdp = m.solve(gap_tolerance=1e-4)
print("status =", res_gdp.status)
print("F =", round(float(res_gdp.objective), 6))
print("x, y, mu =", round(value(res_gdp, "x"), 6),
round(value(res_gdp, "y"), 6), round(value(res_gdp, "mu"), 6))
# The certified global optimum: x=2, y=2, F=2, with the cap y<=x active (mu=2).
assert res_gdp.status == "optimal"
assert abs(res_gdp.objective - 2.0) < 1e-4
assert abs(value(res_gdp, "x") - 2.0) < 1e-4 and abs(value(res_gdp, "y") - 2.0) < 1e-4
nodes_gdp = res_gdp.node_count
print("nodes =", nodes_gdp)
status = optimal
F = 2.0
x, y, mu = 2.0 2.0 2.0
nodes = 3
Contrast: the smooth bilinear encoding#
The historical alternative is to keep complementary slackness as the smooth bilinear equality \(\mu \cdot g = 0\) and let the spatial solver handle the nonconvex product directly. It reaches the same certified optimum, but its McCormick relaxation of \(\mu\,g\) is loose — the solver must spatially branch to tighten it — whereas the disjunction branches on a single binary.
On a problem this small the difference is a couple of nodes, not a dramatic
speedup; the point is the mechanism, not the margin. (Measured, the margin here
is also flat in the multiplier bound: sweeping mu’s upper bound from 20 to
10,000 leaves both encodings at 3 and 5 nodes respectively, so the advantage on
this model is structural rather than something that grows as the box loosens.)
m, x, y, mu, g = build_kkt_model()
m.subject_to(mu * (-g) == 0, name="compl_slack") # smooth bilinear mu * g == 0
# The bilinear MILP relaxation drives the simplex ratio test through near-zero
# pivots, so NumPy flags a benign overflow (the resulting `inf` ratios are
# discarded downstream). Silence it so the teaching output stays clean.
with np.errstate(over="ignore", divide="ignore", invalid="ignore"):
res_bil = m.solve(gap_tolerance=1e-4)
print("status =", res_bil.status)
print("F =", round(float(res_bil.objective), 6))
print("x, y, mu =", round(value(res_bil, "x"), 6),
round(value(res_bil, "y"), 6), round(value(res_bil, "mu"), 6))
assert res_bil.status == "optimal"
assert abs(res_bil.objective - 2.0) < 1e-4 # same certified optimum
nodes_bil = res_bil.node_count
print()
print(f"nodes (gdp disjunction) = {nodes_gdp}")
print(f"nodes (smooth bilinear) = {nodes_bil}")
print(f"the disjunction explores {nodes_bil - nodes_gdp} fewer nodes "
f"({nodes_bil / max(nodes_gdp, 1):.2f}x)")
INFO pounce_algorithm::conv_check::opt_error: certifying with a dual infeasibility above dual_inf_tol: it is within the scale-relative floor set by the terms the Lagrangian gradient is built from (dual_inf_scale_kappa=0 disables) dual_inf=22.837890625 dual_scale=9967700070716.86 dual_inf_tol=1.0 bound=99677.0007071686
status = optimal
F = 2.0
x, y, mu = 2.0 2.0 2.0
nodes (gdp disjunction) = 3
nodes (smooth bilinear) = 5
the disjunction explores 2 fewer nodes (1.67x)
Takeaways#
The KKT reduction turns the bilevel program into a single-level MPEC; its equilibrium content lives entirely in the complementary-slackness pairs.
Model.complementarity()encodes each pair as the exact disjunction \([\mu = 0] \vee [g = 0]\), so branch-and-bound branches the finite either/or and the integrality-aware bound tightening infers one side is zero when the other is bounded away from it — the disjunctive advantage over the smooth \(\mu\,g = 0\) product [Fortuny-Amat and McCarl, 1981, Raman and Grossmann, 1994].Both encodings certify the same global optimum \(F^\star = 2\); the disjunction gets there in fewer nodes (3 vs 5 here). The margin is small because the model is tiny — what matters is that the disjunction branches a finite either/or instead of spatially branching a bilinear product.
For larger models, discopt.bilevel.BilevelProblem automates exactly this
pipeline — building the KKT system, folding follower bounds, and lowering the
complementary-slackness pairs through the same complementarity() machinery.
The manual build_kkt_model() above is the transparent version of what it does
under the hood; see the bilevel notebook for the automated
front-end and its convexity gate.
References#
See the References page for full citations.