Bilevel Optimization (Stackelberg games)

Bilevel Optimization (Stackelberg games)#

A bilevel program models a leader who optimizes subject to a follower solving its own optimization problem in response:

\[\min_{x,y}\; F(x,y)\quad\text{s.t.}\quad G(x,y)\le 0,\quad y \in \arg\min_{y'}\{ f(x,y') : g(x,y') \le 0 \}.\]

The leader picks \(x\); the follower responds with the \(y\) that is optimal for its problem. This models toll setting, network interdiction, pricing, and Stackelberg competition [Bard, 1998, Colson et al., 2007, Dempe, 2002].

The KKT reduction#

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 we can replace y ∈ argmin{…} by

stationarity \(\nabla_y L = 0\) with \(L = f + \sum_i \mu_i g_i\),
primal feasibility \(g_i \le 0\),
dual feasibility \(\mu_i \ge 0\),
complementarity \(0 \le \mu_i \perp -g_i \ge 0\).

The result is a single-level MPEC, which discopt already solves (discopt.mpec: GDP / SOS1 / Scholtes). discopt’s BilevelProblem builds this reduction automatically. It handles the optimistic case; the lower level must be convex in \(y\) — an LP, a convex-QP, or a convex-NLP proved convex by an interval-Hessian test over the box (so the natural linearly-coupled form such as \(e^{y} - x\,y\) is accepted). Nonconvex / integer-follower / pessimistic problems are refused loudly rather than silently mis-solved [Kleinert et al., 2021].

Stationarity needs a symbolic derivative#

So that \(\nabla_y L = 0\) stays inside discopt’s certified global path, the gradient must be an ordinary expression (not an opaque numeric callback). The discopt.bilevel.symbolic_diff engine provides diff/grad over the expression DAG:

import discopt.modeling as dm
from discopt.bilevel import diff

m = dm.Model('demo')
x = m.continuous('x', lb=-5, ub=5)
expr = x**2 + 3 * x
print('d/dx (x^2 + 3x) =', diff(expr, x))   # simplifies to (x + x) + 3

d/dx (x^2 + 3x) = ((2 * x) + 3)

Worked example: toll setting#

A road authority (leader) sets a toll \(t\) on an arc to maximize revenue \(t\,f\). Drivers (follower) choose the flow \(f\) to minimize a convex congestion-plus-toll cost \(\tfrac12 f^2 - 5f + t f\). The follower is a convex QP in \(f\), so BilevelProblem replaces it by its KKT conditions. We build the follower’s KKT system with build_kkt_system() and inspect the emitted conditions (pure Python — no solver needed for this step):

from discopt.bilevel import BilevelProblem

m = dm.Model('toll')
t = m.continuous('toll', lb=0, ub=5)   # leader: toll
f = m.continuous('flow', lb=0, ub=8)   # follower: flow
m.minimize(-(t * f))                    # maximize revenue t*f

bl = BilevelProblem(
    m, upper_vars=[t], lower_vars=[f],
    lower_objective=0.5 * f * f - 5 * f + t * f,  # convex QP in f
    lower_constraints=[],                          # only the flow bounds (folded in)
)

# build_kkt_system() emits the follower's KKT conditions (sound math) and returns
# them for inspection, independent of how the complementarity is later encoded.
bl.build_kkt_system()

print('stationarity:', bl.kkt.stationarity[0].body)
for p in bl.kkt.comp_pairs:
    print('complementarity: 0 <=', p.f, 'perp', p.g, '>= 0')

stationarity: ((((((0.5 * flow) + (0.5 * flow)) - 5) + toll) + (bl_mu0 * -1)) + bl_mu1)
complementarity: 0 <= bl_mu0 perp neg((0 - flow)) >= 0
complementarity: 0 <= bl_mu1 perp neg((flow - 8)) >= 0

The follower’s finite bounds (\(0 \le f \le 8\)) are folded into the KKT system automatically — they are genuine follower constraints, so omitting them would be wrong whenever the follower optimum sits on a bound.

Solving and the optimum#

Strong duality. Call bl.formulate(method='strong_duality') then m.solve(). Strong duality keeps stationarity, primal and dual feasibility, and replaces the per-row complementarity with the single aggregate equality \(\sum_i \mu_i g_i = 0\) (one bilinear equality, no disjunction).

We can check the answer by hand: the follower’s optimum is \(f = 5 - t\), so the leader’s revenue is \(t(5-t)\), maximized at \(t^\* = 2.5\), \(f^\* = 2.5\), revenue \(6.25\).

That equality is nonconvex, so certification is not guaranteed — but nor is it precluded, since the spatial branch-and-bound can still close the gap. On a model this small it does, and the cell below reports what actually came back rather than predicting it.

Big-M KKT encodings. mpec_method='gdp'/'sos1' need a finite bound on each follower multiplier: a KKT multiplier has no a-priori upper bound, and a sentinel-sized big-M is numerically vacuous — it would let the solver certify a follower-infeasible point (a false optimum), so BilevelProblem refuses the big-M methods loudly for unbounded multipliers. Deriving a valid bound automatically is NP-hard in general [Kleinert et al., 2021], so if you can assert one for your problem, pass it explicitly as multiplier_ub.

The bound is a modeling assertion, exactly like a user-supplied big-M: it must be \(\ge\) the largest multiplier at the true follower optimum for every leader decision, or it silently cuts the true response. solve() warns if a multiplier ends up at its supplied bound.

def build_toll(multiplier_ub=None):
    """Rebuild the toll model; pass multiplier_ub to enable the big-M encodings."""
    mm = dm.Model("toll")
    tt = mm.continuous("toll", lb=0, ub=5)
    ff = mm.continuous("flow", lb=0, ub=8)
    mm.minimize(-(tt * ff))
    kwargs = {} if multiplier_ub is None else {"multiplier_ub": multiplier_ub}
    prob = BilevelProblem(
        mm, upper_vars=[tt], lower_vars=[ff],
        lower_objective=0.5 * ff * ff - 5 * ff + tt * ff,
        lower_constraints=[], **kwargs,
    )
    return mm, prob


ANALYTIC_REVENUE = 6.25  # t*(5-t) at t = 2.5

# 1. Strong duality (no multiplier bound needed).
mm, prob = build_toll()
prob.formulate(method="strong_duality")
r = mm.solve()
print(f"strong_duality      : toll={r.x['toll']:.4f}  flow={r.x['flow']:.4f}  "
      f"revenue={-r.objective:.4f}  gap_certified={r.gap_certified}")
assert abs(-r.objective - ANALYTIC_REVENUE) < 1e-6, -r.objective

# 2. Big-M encodings WITHOUT a multiplier bound: these must refuse, not guess.
mm, prob = build_toll()
try:
    prob.formulate(method="kkt", mpec_method="gdp")
    print("gdp (no bound)      : UNEXPECTED -- should have refused")
except NotImplementedError as exc:
    print(f"gdp (no bound)      : refused, as designed -- {str(exc).split(chr(10))[0][:70]}...")

# 3. Big-M encodings WITH an asserted bound.
for mpec_method in ("gdp", "sos1"):
    mm, prob = build_toll(multiplier_ub=50.0)
    prob.formulate(method="kkt", mpec_method=mpec_method)
    r = mm.solve()
    print(f"kkt/{mpec_method:<4s} (ub=50)    : toll={r.x['toll']:.4f}  flow={r.x['flow']:.4f}  "
          f"revenue={-r.objective:.4f}  gap_certified={r.gap_certified}")
    assert abs(-r.objective - ANALYTIC_REVENUE) < 1e-6, -r.objective

print(f"\nAll solved arms match the analytic optimum ({ANALYTIC_REVENUE}).")

/Users/jkitchin/projects/discopt/python/discopt/debug/__init__.py:98: UserWarning: Variables with very large or infinite declared bounds: bl_mu0 (lb=0, ub=9.999e+19 [default]), bl_mu1 (lb=0, ub=9.999e+19 [default]). A bound marked [default] is the box applied to a column you declared with no bounds: 9.999e+19, which is FINITE (it sits just below the 1e+20 infinity sentinel), so the solve can return a certified 'optimal' sitting on that corner rather than 'unbounded'. NLP solvers may fail (NaN, iteration_limit) when bounds exceed ~1e15. Add tighter explicit bounds, e.g. m.continuous('x', lb=0, ub=1000).
  return fn(*args, **kwargs)
 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=1.4667778015136719 dual_scale=15684713060.164063 dual_inf_tol=1.0 bound=156.84713060164063
 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=2.850830078125 dual_scale=1771142500419.375 dual_inf_tol=1.0 bound=17711.425004193752
 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=40.04273986816406 dual_scale=110702145987.36311 dual_inf_tol=1.0 bound=1107.021459873631
 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=8.390632629394531 dual_scale=18621633243.747612 dual_inf_tol=1.0 bound=186.21633243747613

strong_duality      : toll=2.5000  flow=2.5000  revenue=6.2500  gap_certified=True
gdp (no bound)      : refused, as designed -- the KKT complementarity multiplier(s) ['bl_mu0', 'bl_mu1'] are unbound...
kkt/gdp  (ub=50)    : toll=2.5000  flow=2.5000  revenue=6.2500  gap_certified=True
kkt/sos1 (ub=50)    : toll=2.5000  flow=2.5000  revenue=6.2500  gap_certified=True

All solved arms match the analytic optimum (6.25).

References#

See the References page for full citations.