0
# Mathematical Programming with Equilibrium Constraints
1

2
Components for modeling mathematical programs with equilibrium constraints (MPEC). MPEC problems involve optimization problems subject to constraints that are themselves defined by the solution of other optimization problems, commonly arising in game theory, economics, and engineering design.
3

4
## Capabilities
5

6
### Complementarity Constraints
7

8
Components for defining complementarity relationships between variables and constraints in equilibrium problems.
9

10
```python { .api }
11
class Complementarity:
12
    """
13
    Complementarity constraint component for MPEC problems.
14
    
15
    Represents complementarity relationships of the form:
16
    0 ≤ x ⊥ f(x) ≥ 0 (x and f(x) cannot both be positive)
17
    """
18
    def __init__(self, *args, **kwargs): ...
19
    
20
    def activate(self):
21
        """Activate this complementarity constraint."""
22
    
23
    def deactivate(self):
24
        """Deactivate this complementarity constraint."""
25
    
26
    def is_active(self):
27
        """Check if complementarity constraint is active."""
28

29
class ComplementarityList:
30
    """
31
    Complementarity constraint list for dynamic constraint creation.
32
    
33
    Allows adding complementarity constraints dynamically during
34
    model construction or solution process.
35
    """
36
    def __init__(self, **kwargs): ...
37
    
38
    def add(self, expr1, expr2):
39
        """
40
        Add complementarity constraint.
41
        
42
        Args:
43
            expr1: First expression in complementarity
44
            expr2: Second expression in complementarity
45
        """
46
```
47

48
### Complementarity Relationship Functions
49

50
Helper functions for defining complementarity relationships between variables and expressions.
51

52
```python { .api }
53
def complements(expr1, expr2):
54
    """
55
    Define complementarity relationship between two expressions.
56
    
57
    Creates a complementarity constraint of the form:
58
    expr1 ≥ 0, expr2 ≥ 0, expr1 * expr2 = 0
59
    
60
    Args:
61
        expr1: First expression (typically variable or constraint)
62
        expr2: Second expression (typically constraint or variable)
63
        
64
    Returns:
65
        Complementarity: Complementarity constraint object
66
    """
67
```
68

69
## Usage Examples
70

71
### Basic Complementarity Problem
72

73
```python
74
from pyomo.environ import *
75
from pyomo.mpec import Complementarity, complements
76

77
model = ConcreteModel()
78

79
# Variables
80
model.x = Var(bounds=(0, None))  # x ≥ 0
81
model.y = Var(bounds=(0, None))  # y ≥ 0
82

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

86
# Complementarity constraint: x ⊥ (x + y - 1)
87
# Either x = 0 or (x + y - 1) = 0 (or both)
88
model.complementarity = Complementarity(
89
    expr=complements(model.x, model.x + model.y - 1)
90
)
91

92
# Alternative syntax
93
model.comp_alt = Complementarity(
94
    expr=complements(model.x >= 0, model.x + model.y - 1 >= 0)
95
)
96
```
97

98
### Variational Inequality Problem
99

100
```python
101
from pyomo.environ import *
102
from pyomo.mpec import Complementarity, complements
103

104
model = ConcreteModel()
105

106
# Decision variables
107
model.x1 = Var(bounds=(0, 10))
108
model.x2 = Var(bounds=(0, 10))
109

110
# Lagrange multipliers for bound constraints
111
model.lambda1 = Var(bounds=(0, None))
112
model.lambda2 = Var(bounds=(0, None))
113
model.mu1 = Var(bounds=(0, None))
114
model.mu2 = Var(bounds=(0, None))
115

116
# KKT conditions as complementarity constraints
117
# Stationarity: ∇f(x) - λ + μ = 0
118
model.stationarity1 = Constraint(
119
    expr=2*model.x1 - model.lambda1 + model.mu1 == 0
120
)
121
model.stationarity2 = Constraint(
122
    expr=2*model.x2 - model.lambda2 + model.mu2 == 0
123
)
124

125
# Complementarity for lower bounds: x ⊥ λ
126
model.comp_lower1 = Complementarity(
127
    expr=complements(model.x1, model.lambda1)
128
)
129
model.comp_lower2 = Complementarity(
130
    expr=complements(model.x2, model.lambda2)
131
)
132

133
# Complementarity for upper bounds: (10 - x) ⊥ μ
134
model.comp_upper1 = Complementarity(
135
    expr=complements(10 - model.x1, model.mu1)
136
)
137
model.comp_upper2 = Complementarity(
138
    expr=complements(10 - model.x2, model.mu2)
139
)
140

141
# Dummy objective (feasibility problem)
142
model.obj = Objective(expr=0)
143
```
144

145
### Bilevel Optimization as MPEC
146

147
```python
148
from pyomo.environ import *
149
from pyomo.mpec import Complementarity, complements
150

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

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

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

160
# Lower-level Lagrange multipliers
161
model.lambda_lower = Var(bounds=(0, None))
162
model.mu_y_lower = Var(bounds=(0, None))
163
model.mu_y_upper = Var(bounds=(0, None))
164

165
# Upper-level objective
166
model.upper_obj = Objective(
167
    expr=(model.x - 3)**2 + (model.y - 2)**2, 
168
    sense=minimize
169
)
170

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

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

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

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

188
# Complementarity for y bounds: y ⊥ μ_y_lower, (5-y) ⊥ μ_y_upper
189
model.comp_y_lower = Complementarity(
190
    expr=complements(model.y, model.mu_y_lower)
191
)
192
model.comp_y_upper = Complementarity(
193
    expr=complements(5 - model.y, model.mu_y_upper)
194
)
195
```
196

197
### Economic Equilibrium Model
198

199
```python
200
from pyomo.environ import *
201
from pyomo.mpec import Complementarity, complements
202

203
model = ConcreteModel()
204

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

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

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

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

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

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

225
# Complementarity: profit maximization conditions
226
# Either q[i] = 0 or price equals marginal cost
227
def complementarity_rule(model, i):
228
    return complements(
229
        model.q[i], 
230
        model.p[i] - model.marginal_cost[i]
231
    )
232

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

236
# Objective (welfare maximization)
237
model.welfare = Objective(
238
    expr=sum(
239
        model.demand_intercept[i] * model.q[i] - 0.5 * model.demand_slope[i] * model.q[i]**2 
240
        - model.marginal_cost[i] * model.q[i]
241
        for i in model.MARKETS
242
    ),
243
    sense=maximize
244
)
245
```
246

247
### Solving MPEC Problems
248

249
```python
250
from pyomo.environ import *
251
from pyomo.mpec import Complementarity, complements
252

253
# Create MPEC model (using previous examples)
254
model = create_mpec_model()
255

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

259
# Smooth complementarity reformulation
260
transform = TransformationFactory('mpec.simple_nonlinear')
261
transform.apply_to(model)
262

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

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

273
# Access solution
274
if results.solver.termination_condition == TerminationCondition.optimal:
275
    print("Solution found:")
276
    for var in model.component_objects(Var):
277
        if var.is_indexed():
278
            for index in var.index():
279
                print(f"{var.name}[{index}] = {value(var[index])}")
280
        else:
281
            print(f"{var.name} = {value(var)}")
282
```
283

284
### Multi-Leader-Follower Game
285

286
```python
287
from pyomo.environ import *
288
from pyomo.mpec import Complementarity, complements
289

290
model = ConcreteModel()
291

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

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

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

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

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

310
model.leader_obj = Objective(
311
    expr=sum(leader_obj_rule(model, i) for i in model.LEADERS),
312
    sense=minimize
313
)
314

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

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

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

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

328
# Follower complementarity conditions
329
def follower_comp_constraint_rule(model, j):
330
    return complements(
331
        model.lambda_f[j], 
332
        6 - sum(model.x[i] for i in model.LEADERS) - model.y[j]
333
    )
334

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

338
model.follower_comp_constraint = Complementarity(model.FOLLOWERS, rule=follower_comp_constraint_rule)
339
model.follower_comp_bound = Complementarity(model.FOLLOWERS, rule=follower_comp_bound_rule)
340
```