Mathematical Programming with Equilibrium Constraints

class Complementarity:
    """
    Complementarity constraint component for MPEC problems.
    
    Represents complementarity relationships of the form:
    0 ≤ x ⊥ f(x) ≥ 0 (x and f(x) cannot both be positive)
    """
    def __init__(self, *args, **kwargs): ...
    
    def activate(self):
        """Activate this complementarity constraint."""
    
    def deactivate(self):
        """Deactivate this complementarity constraint."""
    
    def is_active(self):
        """Check if complementarity constraint is active."""

class ComplementarityList:
    """
    Complementarity constraint list for dynamic constraint creation.
    
    Allows adding complementarity constraints dynamically during
    model construction or solution process.
    """
    def __init__(self, **kwargs): ...
    
    def add(self, expr1, expr2):
        """
        Add complementarity constraint.
        
        Args:
            expr1: First expression in complementarity
            expr2: Second expression in complementarity
        """

def complements(expr1, expr2):
    """
    Define complementarity relationship between two expressions.
    
    Creates a complementarity constraint of the form:
    expr1 ≥ 0, expr2 ≥ 0, expr1 * expr2 = 0
    
    Args:
        expr1: First expression (typically variable or constraint)
        expr2: Second expression (typically constraint or variable)
        
    Returns:
        Complementarity: Complementarity constraint object
    """

from pyomo.environ import *
from pyomo.mpec import Complementarity, complements

model = ConcreteModel()

# Variables
model.x = Var(bounds=(0, None))  # x ≥ 0
model.y = Var(bounds=(0, None))  # y ≥ 0

# Objective
model.obj = Objective(expr=model.x + model.y, sense=minimize)

# Complementarity constraint: x ⊥ (x + y - 1)
# Either x = 0 or (x + y - 1) = 0 (or both)
model.complementarity = Complementarity(
    expr=complements(model.x, model.x + model.y - 1)
)

# Alternative syntax
model.comp_alt = Complementarity(
    expr=complements(model.x >= 0, model.x + model.y - 1 >= 0)
)

from pyomo.environ import *
from pyomo.mpec import Complementarity, complements

model = ConcreteModel()

# Decision variables
model.x1 = Var(bounds=(0, 10))
model.x2 = Var(bounds=(0, 10))

# Lagrange multipliers for bound constraints
model.lambda1 = Var(bounds=(0, None))
model.lambda2 = Var(bounds=(0, None))
model.mu1 = Var(bounds=(0, None))
model.mu2 = Var(bounds=(0, None))

# KKT conditions as complementarity constraints
# Stationarity: ∇f(x) - λ + μ = 0
model.stationarity1 = Constraint(
    expr=2*model.x1 - model.lambda1 + model.mu1 == 0
)
model.stationarity2 = Constraint(
    expr=2*model.x2 - model.lambda2 + model.mu2 == 0
)

# Complementarity for lower bounds: x ⊥ λ
model.comp_lower1 = Complementarity(
    expr=complements(model.x1, model.lambda1)
)
model.comp_lower2 = Complementarity(
    expr=complements(model.x2, model.lambda2)
)

# Complementarity for upper bounds: (10 - x) ⊥ μ
model.comp_upper1 = Complementarity(
    expr=complements(10 - model.x1, model.mu1)
)
model.comp_upper2 = Complementarity(
    expr=complements(10 - model.x2, model.mu2)
)

# Dummy objective (feasibility problem)
model.obj = Objective(expr=0)

from pyomo.environ import *
from pyomo.mpec import Complementarity, complements

# Upper-level and lower-level variables
model = ConcreteModel()

# Upper-level decision variables
model.x = Var(bounds=(0, 10))

# Lower-level decision variables  
model.y = Var(bounds=(0, 5))

# Lower-level Lagrange multipliers
model.lambda_lower = Var(bounds=(0, None))
model.mu_y_lower = Var(bounds=(0, None))
model.mu_y_upper = Var(bounds=(0, None))

# Upper-level objective
model.upper_obj = Objective(
    expr=(model.x - 3)**2 + (model.y - 2)**2, 
    sense=minimize
)

# Lower-level KKT conditions
# Stationarity: ∇_y f_lower(x,y) - λ∇_y g(x,y) + μ_lower - μ_upper = 0
model.ll_stationarity = Constraint(
    expr=2*(model.y - 1) - model.lambda_lower * 1 + model.mu_y_lower - model.mu_y_upper == 0
)

# Lower-level constraint: x + y <= 4
model.ll_constraint = Constraint(expr=model.x + model.y <= 4)

# Complementarity conditions for lower-level problem
# Primal feasibility and dual feasibility handled by bounds and constraints

# Complementary slackness: λ ⊥ (4 - x - y) 
model.comp_ll_constraint = Complementarity(
    expr=complements(model.lambda_lower, 4 - model.x - model.y)
)

# Complementarity for y bounds: y ⊥ μ_y_lower, (5-y) ⊥ μ_y_upper
model.comp_y_lower = Complementarity(
    expr=complements(model.y, model.mu_y_lower)
)
model.comp_y_upper = Complementarity(
    expr=complements(5 - model.y, model.mu_y_upper)
)

from pyomo.environ import *
from pyomo.mpec import Complementarity, complements

model = ConcreteModel()

# Markets
model.MARKETS = Set(initialize=['A', 'B', 'C'])

# Production variables (quantity produced in each market)
model.q = Var(model.MARKETS, bounds=(0, 100))

# Price variables
model.p = Var(model.MARKETS, bounds=(0, None))

# Cost parameters
model.marginal_cost = Param(model.MARKETS, initialize={'A': 2, 'B': 3, 'C': 4})

# Demand parameters  
model.demand_intercept = Param(model.MARKETS, initialize={'A': 50, 'B': 40, 'C': 60})
model.demand_slope = Param(model.MARKETS, initialize={'A': 1, 'B': 0.8, 'C': 1.2})

# Market clearing: supply equals demand
def market_clearing_rule(model, i):
    return model.q[i] == max(0, model.demand_intercept[i] - model.demand_slope[i] * model.p[i])

# Complementarity: profit maximization conditions
# Either q[i] = 0 or price equals marginal cost
def complementarity_rule(model, i):
    return complements(
        model.q[i], 
        model.p[i] - model.marginal_cost[i]
    )

model.market_clearing = Constraint(model.MARKETS, rule=market_clearing_rule)
model.profit_max = Complementarity(model.MARKETS, rule=complementarity_rule)

# Objective (welfare maximization)
model.welfare = Objective(
    expr=sum(
        model.demand_intercept[i] * model.q[i] - 0.5 * model.demand_slope[i] * model.q[i]**2 
        - model.marginal_cost[i] * model.q[i]
        for i in model.MARKETS
    ),
    sense=maximize
)

from pyomo.environ import *
from pyomo.mpec import Complementarity, complements

# Create MPEC model (using previous examples)
model = create_mpec_model()

# Method 1: Transform to smooth NLP using complementarity reformulation
from pyomo.mpec import TransformationFactory

# Smooth complementarity reformulation
transform = TransformationFactory('mpec.simple_nonlinear')
transform.apply_to(model)

solver = SolverFactory('ipopt')
results = solver.solve(model, tee=True)

# Method 2: Use specialized MPEC solver (if available)
try:
    mpec_solver = SolverFactory('knitro_mpec')  # Commercial solver
    results = mpec_solver.solve(model, tee=True)
except:
    print("MPEC solver not available, using NLP reformulation")

# Access solution
if results.solver.termination_condition == TerminationCondition.optimal:
    print("Solution found:")
    for var in model.component_objects(Var):
        if var.is_indexed():
            for index in var.index():
                print(f"{var.name}[{index}] = {value(var[index])}")
        else:
            print(f"{var.name} = {value(var)}")

from pyomo.environ import *
from pyomo.mpec import Complementarity, complements

model = ConcreteModel()

# Players
model.LEADERS = Set(initialize=[1, 2])
model.FOLLOWERS = Set(initialize=[1, 2, 3])

# Leader decision variables
model.x = Var(model.LEADERS, bounds=(0, 10))

# Follower decision variables
model.y = Var(model.FOLLOWERS, bounds=(0, 5))

# Follower multipliers for their KKT conditions
model.lambda_f = Var(model.FOLLOWERS, bounds=(0, None))
model.mu_f = Var(model.FOLLOWERS, bounds=(0, None))

# Leader objectives (Stackelberg competition)
def leader_obj_rule(model, i):
    return (model.x[i] - 2)**2 + sum(model.y[j] for j in model.FOLLOWERS) / len(model.FOLLOWERS)

model.leader_obj = Objective(
    expr=sum(leader_obj_rule(model, i) for i in model.LEADERS),
    sense=minimize
)

# Follower KKT conditions
def follower_stationarity_rule(model, j):
    # Each follower solves: min (y_j - 1)^2 subject to sum(x_i) + y_j <= 6
    return 2*(model.y[j] - 1) + model.lambda_f[j] - model.mu_f[j] == 0

model.follower_stationarity = Constraint(model.FOLLOWERS, rule=follower_stationarity_rule)

# Follower feasibility constraint
def follower_constraint_rule(model, j):
    return sum(model.x[i] for i in model.LEADERS) + model.y[j] <= 6

model.follower_constraints = Constraint(model.FOLLOWERS, rule=follower_constraint_rule)

# Follower complementarity conditions
def follower_comp_constraint_rule(model, j):
    return complements(
        model.lambda_f[j], 
        6 - sum(model.x[i] for i in model.LEADERS) - model.y[j]
    )

def follower_comp_bound_rule(model, j):
    return complements(model.y[j], model.mu_f[j])

model.follower_comp_constraint = Complementarity(model.FOLLOWERS, rule=follower_comp_constraint_rule)
model.follower_comp_bound = Complementarity(model.FOLLOWERS, rule=follower_comp_bound_rule)