Complementarity and MPEC Models

Complementarity and MPEC Models#

Many equilibrium, bilevel, and contact problems include complementarity conditions

\[ 0 \le f(x) \;\perp\; g(x) \ge 0 \qquad\Longleftrightarrow\qquad f \ge 0,\; g \ge 0,\; f\cdot g = 0. \]

A program with such constraints is a Mathematical Program with Equilibrium Constraints (MPEC). The product condition \(f\cdot g = 0\) is nonconvex and violates the standard constraint qualifications, so a naive NLP solve stalls [Luo et al., 1996].

The fluent front-end is Model.complementarity(x, y), and discopt.mpec provides three reformulations behind it:

GDP disjunction (the default) — encode the either/or structure exactly as a generalized disjunction \([x = 0] \vee [y = 0]\) and let discopt’s GDP machinery reformulate it (big-M / hull) into the MILP/MINLP relaxation [Raman and Grossmann, 1994]. This is exact and pairs naturally with the spatial branch-and-bound.
Scholtes regularization [Scholtes, 2001] — relax \(f\cdot g = 0\) to \(f\cdot g \le t\) and drive \(t \to 0\) through a homotopy of local NLP solves; fast for smooth continuous MPECs.
Exact SOS1 — place \(\{x, y\}\) in a special-ordered-set-1 so at most one is nonzero; the global solver branches on the two disjuncts.

All three are reachable both through the fluent Model.complementarity(x, y, method=...) and through solve_mpec(model, [...], method=...).

import os
os.environ["JAX_PLATFORMS"] = "cpu"
os.environ["JAX_ENABLE_X64"] = "1"

import numpy as np
import discopt.modeling.core as dm
from discopt.mpec import complementarity, solve_mpec, tighten_complementarity_bounds

A reference MPEC#

We minimize the distance from \((x, y)\) to \((1, 1)\) subject to \(0 \le x \perp y \ge 0\). The complementarity forces \(x\cdot y = 0\), so the global optima sit at \((1, 0)\) and \((0, 1)\), both with objective \(1\).

def distance_model():
    m = dm.Model("mpec_distance")
    x = m.continuous("x", lb=0, ub=10)
    y = m.continuous("y", lb=0, ub=10)
    m.minimize((x - 1) ** 2 + (y - 1) ** 2)
    return m, x, y

The fluent API and the default GDP reformulation#

Model.complementarity(x, y) records the condition on the model and reformulates it immediately, so a plain model.solve() certifies the result. With no method given it uses method="gdp" — the disjunction \([x=0]\vee[y=0]\) — which is exact and solved to global optimality by branch-and-bound.

m, x, y = distance_model()
m.complementarity(x, y)          # method="gdp" by default
res_gdp = m.solve()
print("status    =", res_gdp.status)
print("objective =", round(float(res_gdp.objective), 6))
print("x, y =", float(np.asarray(res_gdp.x["x"])), float(np.asarray(res_gdp.x["y"])))

status    = optimal
objective = 1.0
x, y = 9.980413130885437e-08 1.0000000023418656

The disjunction is satisfied at \((0, 1)\) with objective \(1\) — a certified global optimum. The same reformulation is available explicitly as solve_mpec(m, [complementarity(x, y)], method="gdp") or discopt.mpec.reformulate_gdp when you want to inspect the generated disjunction.

Scholtes regularization#

The homotopy solves a sequence of smooth NLPs with \(t = 1, 10^{-1}, \dots\), warm-starting each from the previous solution. As \(t \to 0\) the iterate approaches the complementary set.

m, x, y = distance_model()
res = solve_mpec(m, [complementarity(x, y)], method="scholtes")
print("x, y =", np.round(res.x, 6))
print("objective =", round(float(res.objective), 6))
print("x*y =", float(res.x[0] * res.x[1]))

x, y = [1. 0.]
objective = 1.0
x*y = 1.66549889400699e-08

The iterate converges to \((0, 1)\) (up to the regularization tolerance) with objective \(1\) and \(x\cdot y \approx 0\) — a global optimum of the MPEC.

Exact SOS1 reformulation#

For a certifiable answer we encode the complementarity exactly. The SOS1 set \(\{x, y\}\) allows at most one member to be nonzero; the global MINLP solver branches on the two disjuncts \(x = 0\) and \(y = 0\).

m2, x2, y2 = distance_model()
res2 = solve_mpec(m2, [complementarity(x2, y2)], method="sos1")
print("status   =", res2.status)
print("objective =", round(float(res2.objective), 6))
# Model.solve returns a name -> value mapping.
print("x, y =", float(np.asarray(res2.x["x"])), float(np.asarray(res2.x["y"])))

status   = optimal
objective = 1.0
x, y = 1.0783785581023457e-07 1.0000000043086208

Both reformulations reach objective \(1\), confirming the same global optimum by two independent routes.

Complementarity bound propagation#

When one side of \(0 \le x \perp y \ge 0\) is provably strictly positive (\(\text{lb} > 0\)), the other side must be zero. This sound implication is exact for single-variable sides and tightens the model before search.

m3 = dm.Model("bp")
a = m3.continuous("a", lb=0.5, ub=3.0)   # strictly positive lower bound
b = m3.continuous("b", lb=0.0, ub=3.0)
n_fixed = tighten_complementarity_bounds(m3, [complementarity(a, b)])
print("variables fixed to zero:", n_fixed)
print("b bounds:", float(b.lb), float(b.ub))

variables fixed to zero: 1
b bounds: 0.0 0.0

When to use which#

Method	Path	Best for
`gdp` (default)	exact disjunction → MILP/MINLP branch-and-bound	exact certificates; the general-purpose choice
`scholtes`	homotopy of local NLPs	smooth continuous MPECs; fast, local
`sos1`	SOS1 global branch-and-bound	exact certificates; small/medium discrete blocks

The GDP disjunction is the default because it is exact and reuses discopt’s branch-and-bound directly. Scholtes is the fast local workhorse for large continuous equilibrium models; SOS1 is an alternative exact encoding that pairs naturally with branch-and-bound when the complementary blocks are modest in number.