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:

  1. 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.

  2. 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.

  3. 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 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 = 1.0000000016600155 6.064465117354872e-11

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 — a local result#

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.

Every method returns the same type — a SolveResult whose x is a {name: value} mapping — but they do not make the same claim. The GDP and SOS1 arms encode the condition exactly and return a certified status="optimal". Scholtes relaxes \(f\cdot g = 0\) to \(f\cdot g \le t\) and solves that locally, so it returns status="local_optimal" with gap_certified=False and no bound: a local result may seed an incumbent after independent verification, and may never become a dual bound. If the reported point came from a stage the subsolver did not converge at, the status is the weaker "local_limit" — the iterate is still returned, as a warm start and with its residuals, but local_optimal is reserved for a point that really is stationary.

m, x, y = distance_model()
res = solve_mpec(m, [complementarity(x, y)], method="scholtes")

# `x` is a {name: value} mapping on every arm, so one branch reads any result.
xv = float(np.asarray(res.x["x"]))
yv = float(np.asarray(res.x["y"]))
print("status        =", res.status)
print("gap_certified =", res.gap_certified, "  bound =", res.bound)
print("x, y          =", round(xv, 6), round(yv, 6))
print("objective     =", round(float(res.objective), 6))
print("x*y           =", xv * yv)

status        = local_optimal
gap_certified = False   bound = None
x, y          = 1.0 0.0
objective     = 1.0
x*y           = 1.6654991185456815e-08

The iterate converges to \((0, 1)\) (up to the regularization tolerance) with objective \(1\) and \(x\cdot y \approx 0\). Note what the result does not say: the status is local_optimal, so this is a stationary point of the regularized problem and carries no global claim, even though it happens to be the MPEC’s global optimum here.

Reading the source residuals#

The result carries mpec_report: residuals measured on the source operands the relation declares, separately from the residual of the rows the lowering generated. The two are different numbers and the difference matters — a regularized point can satisfy every generated row while violating the condition the user wrote. Every residual records the formula that produced it, and the trace reports the worst-case residual a Scholtes regularization at \(t\) admits [Scholtes, 2001], so a small residual is never read as exact orthogonality.

That scale is an upper bound on the admitted residual, not a limit on attainable accuracy: \(f\cdot g \le t\) also admits \((f, g) = (0, 1)\), whose residual is exactly \(0\) at every positive \(t\), so landing well below the scale is ordinary. It also depends on which residual definition is read — \(\sqrt{t}\) for min, \(t\) for the product form, \((2-\sqrt{2})\sqrt{t}\) for Fischer–Burmeister — which is why each Residual carries its own admitted_scale beside the formula that produced it.

within_admitted_scale compares the residual against that bound and nothing else — it does not say the relaxed rows hold at the point. At \(t=10^{-8}\), \((f,g)=(10^{-4},1)\) sets the flag while \(f\cdot g\) exceeds \(t\) by four orders of magnitude. Feasibility of the generated rows is the separate lowered_row_residual.

rep = res.mpec_report
print("relations measured:", rep.n_scalar_relations)
print(f"source complementarity : {rep.complementarity.value:.3e}")
print(f"lowered-row residual   : {rep.lowered_row_residual.value:.3e}")
print(f"source primal feasible : {rep.primal_feasibility.value:.3e}")
print("definition:", rep.complementarity.definition)
print()
trace = rep.continuation
print("continuation:", len(trace.stages), "stages,",
      f"final t={trace.final_t:.1e}, {trace.termination_reason},",
      f"converged={trace.converged}")
print(f"admitted min-residual  : {trace.admitted_residual_scale:.1e}",
      "(an upper bound on what the relaxation allows, not a floor)")
print("admitted scale here   :", rep.complementarity.admitted_scale,
      "| within it:", rep.complementarity.within_admitted_scale)
print("lowered-row residual is the separate feasibility number:",
      f"{rep.lowered_row_residual.value:.3e}")

relations measured: 1
source complementarity : 1.665e-08
lowered-row residual   : 6.655e-09
source primal feasible : 0.000e+00
definition: max over scalar relations of r_i = min(f_i(x), g_i(x)); reported as max_i r_i; computed on the SOURCE operands (#1147 provenance), not on the lowered rows

continuation: 10 stages, final t=1.0e-08, t_min_reached, converged=True
admitted min-residual  : 1.0e-04 (an upper bound on what the relaxation allows, not a floor)
admitted scale here   : 0.0001 | within it: True
lowered-row residual is the separate feasibility number: 6.655e-09

On this well-behaved instance the two residuals agree. They come apart exactly where it matters — when the relaxation’s optimum sits away from the complementary set, the generated rows are satisfied while the source condition is not:

# max x + y  s.t.  x == y,  0 <= x _|_ y >= 0.  Complementarity forces x = y = 0,
# but the regularization only asks x*y <= t, so the local optimum is x = y ~ sqrt(t).
mb = dm.Model("balanced")
xb = mb.continuous("x", lb=0, ub=1)
yb = mb.continuous("y", lb=0, ub=1)
mb.subject_to(xb == yb)
mb.minimize(-(xb + yb))

rb = solve_mpec(mb, [complementarity(xb, yb)], method="scholtes", t_min=1e-8)
pb = rb.mpec_report
print(f"objective              : {rb.objective:.3e}   (true MPEC optimum is 0)")
print(f"lowered-row residual   : {pb.lowered_row_residual.value:.3e}   <- rows satisfied")
print(f"source complementarity : {pb.complementarity.value:.3e}   <- condition violated")
print("source_satisfied       :", pb.source_satisfied)

# Which is why the point is refused as an incumbent, and why a local result may
# never become a dual bound: its objective is BETTER than the true optimum.
from discopt.mpec_report import accept_local_incumbent
print("accept_local_incumbent :", accept_local_incumbent(mb, rb))

objective              : -2.828e-04   (true MPEC optimum is 0)
lowered-row residual   : 1.000e-08   <- rows satisfied
source complementarity : 1.414e-04   <- condition violated
source_satisfied       : False
accept_local_incumbent : None

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.0000000015057526 9.264116932692925e-11

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

Vectorized (elementwise) complementarity#

Equilibrium and KKT systems rarely have a single complementarity pair — they have one per constraint or per commodity. complementarity accepts array operands and treats them elementwise [Luo et al., 1996]: for every index \(i\) independently,

\[ 0 \le x_i \;\perp\; y_i \ge 0 \qquad (x_i\,y_i = 0,\; x_i,y_i \ge 0). \]

Each index becomes its own \((x_i=0)\lor(y_i=0)\) disjunction, so the solver may choose a different active side per index. This is strictly different from — and not captured by — a single “all \(x=0\) or all \(y=0\)” choice.

The instance below makes the distinction observable. The objective wants \(y_0=5\) (forcing \(x_0=0\)) for the first pair but \(x_1=5\) (forcing \(y_1=0\)) for the second — a mixed selection. The elementwise optimum is \(0\); an all-or-nothing encoding would be stuck at \(\approx 25\).

n = 2
mv = dm.Model("vector_mpec")
xv = mv.continuous("x", shape=n, lb=0, ub=10)
yv = mv.continuous("y", shape=n, lb=0, ub=10)
# Pair 0 pulls y0 -> 5 (so x0 = 0); pair 1 pulls x1 -> 5 (so y1 = 0).
mv.minimize((yv[0] - 5) ** 2 + (xv[1] - 5) ** 2)
mv.complementarity(xv, yv)          # one call, n elementwise pairs
res = mv.solve()

xs = np.asarray(res.x["x"], dtype=float)
ys = np.asarray(res.x["y"], dtype=float)
print("objective:", round(float(res.objective), 6))
print("x:", np.round(xs, 4))
print("y:", np.round(ys, 4))
print("elementwise products x_i*y_i:", np.round(xs * ys, 6))

objective: 0.0
x: [0. 5.]
y: [5. 0.]
elementwise products x_i*y_i: [0. 0.]

Checking a phase-changing POUNCE flash point#

Issue #1526 uses a small but deliberately difficult physical fixture: POUNCE’s ethane / n-butane flash at 10 bar. Its two source pairs are phase amount against stability slack,

\[0 \le \beta \perp 1-S_y \ge 0, \qquad 0 \le 1-\beta \perp 1-S_x \ge 0.\]

The PR-fast comparison runs 250, 268, 300, 324, and 350 K with gdp, sos1, and scholtes; the benchmark runner can execute POUNCE’s full 34-temperature path. This is a feasibility and regime check, not an objective benchmark: the flash objective is identically zero, so a zero certified gap says only that DiscOpt certified the returned feasible point. The evidence is the full physical point, the regime, and separately defined source residuals.

The exact gdp and sos1 arms are seeded from the POUNCE point: it is passed as initial_solution once accept_local_incumbent accepts it. With a zero objective that incumbent already meets the bound at the root, so the one-node certificates are a feasibility and regime check of POUNCE’s point, not an independent exact cross-validation. Run unseeded, the exact GDP still certifies 300 K and 324 K in one node but reaches a 60 s limit at 250 K after 419 nodes without a certificate; a genuine cold-start comparison is future work.

The global model does not hide Peng–Robinson [Peng and Robinson, 1976] root selection inside a numerical callback. For each phase it declares a SELECT_ONE identity between one-real-root and three-real-root branches. The three-root branch uses all three Vieta identities plus ordering and selects the smallest root for liquid or largest root for vapor. A second explicit disjunction requires at least one \(|\log K_i| \ge 10^{-4}\); without it, the ever-present trivial stationary point \(x=y=z\) can be globally certified because the objective is zero. The comparison artifact records both selected identities.

Why POUNCE is exercised differently here#

The exact Boolean choices belong to GDP/SOS1, not to a local NLP. The Scholtes arm therefore fixes the root-count and nontriviality branches selected by the independently checked POUNCE point, keeps their algebraic equations smooth, and regularizes only the two physical complementarity pairs [Scholtes, 2001] with \(t=10^{-2},10^{-3},\ldots,10^{-8}\). Its result is always local: it may supply a primal incumbent after accept_local_incumbent rechecks the source conditions, but it never supplies a dual bound or certified gap.

This separation also reflects two POUNCE findings from the original fixture. First, an exact-product finish overdetermines this square equilibrium system, so POUNCE’s supported route needs its product-inequality fallback. Second, Cardano’s formula approaches an inverse-trigonometric endpoint near a repeated cubic root and needs forward-over-forward Hessians. DiscOpt’s comparison instead carries the roots as algebraic variables, avoiding that differentiation singularity while exposing the one-root/three-root identity. At an exactly repeated root the two closed discriminant branches overlap; the current benchmark refuses points inside a \(10^{-8}\) discriminant guard rather than pretending it can encode a strict root-count distinction with closed inequalities.

Bubble and dew points are biactive: both sides of one complementarity pair are zero. SELECT_ONE still gives the right feasible union, but either selected branch label is admissible there. Consequently the artifact compares the physical point and records all true active sides; it does not demand a unique branch label at a transition. Source balance, isofugacity, EOS, root-selection, nontriviality, bound, sign, and complementarity residuals remain separate from residuals of the generated lowered rows.

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.