I have a Mixed Complementarity Problem (MCP) that represents a market under perfect competition. The model is written in GAMS and works as expected. But when translating it to Python, using Pyomo, the model is infeasible. Both models use the Path solver.
It appears to be the Pyomo c1_rule constraint that is causing infeasibility. If I change the Pyomo c1_rule relationship to be <= rather than == then the model solves to optimality, but the solution makes no sense.
Any ideas about what I'm doing wrong?
The variables and equations in GAMS are:
positive variables Production(g), Demand; free variables Price; elasticity = -0.3; p0 = 20; g0 = 3000; rho = p0 / elasticity / g0; Production.up(g) = StatData(g, 'GenMax'); equations SRMC(g), PriceCurve, Balance; SRMC(g).. MarginalCost(g) - Price =e= 0; PriceCurve.. Price - (p0 + rho * (Demand - g0)) =e= 0; Balance.. sum(g, Production(g)) - Demand =e= 0; model PerfectComp / SRMC.Production, Balance.Price, PriceCurve.Demand /; solve PerfectComp using MCP;
Equivalent code in Pyomo:
Model.Elasticity = pyo.Param(within = pyo.Reals, initialize = -0.3) Model.pSet = pyo.Param(within = pyo.NonNegativeReals, initialize = 20) Model.qSet = pyo.Param(within = pyo.NonNegativeReals, initialize = 3000) Rho = Model.pSet / Model.Elasticity / Model.qSet Model.Production = pyo.Var(Model.Generators, domain = pyo.NonNegativeReals) Model.Demand = pyo.Var(domain = pyo.NonNegativeReals) Model.Price = pyo.Var(domain = pyo.Reals) def rule_capacity(Model, S): return Model.Production[S] <= Model.GMax[S] Model.MaxCapacity = pyo.Constraint(Model.Generators, rule = rule_capacity) def c1_rule(Model, S): return mpec.complements(Model.VarCost[S] - Model.Price == 0, Model.Production[S] >= 0) def c2_rule(Model): return mpec.complements(Model.Price - (Model.pSet + Rho * (Model.Demand - Model.qSet)) == 0, Model.Demand >= 0) def c3_rule(Model): return mpec.complements(sum(Model.Production[s] for s in Model.Generators) - Model.Demand == 0, Model.Price) Model.c1 = mpec.Complementarity(Model.Generators, rule = c1_rule) Model.c2 = mpec.Complementarity(rule = c2_rule) Model.c3 = mpec.Complementarity(rule = c3_rule) Solver = pyo.SolverFactory('path') Results = Solver.solve(Model, load_solutions = False, tee = True)