3
$\begingroup$

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)
Improve this question
$\endgroup$

5 Answers 5

Reset to default
4
$\begingroup$

I am not sure how an MCP solver/wrapper works, but it seems the solving process, somewhat, is similar to the Lagrangian relaxation method.

In the following picture you can see a simple transportation problem and whose equivalent MCP that implemented in GAMS.

enter image description here

Also, for more details took a look at these links:

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you for your answer. You're right that MCP is similar to Lagrangian relaxation. In deriving constraints, generally we take the first derivative of the terms in a constraint and, in a sense, multiple by the dual variable associated with that constraint. I've done that to create my GAMS model, which works correctly. The problem is that a literal translation to Pyomo results in an infeasible model. I assume that there's something about how Pyomo implements the MCP process that differs from GAMS, so my Pyomo version is incorrectly specified somehow. $\endgroup$
    – Solver Max
    Dec 5, 2022 at 17:51
  • 1
    $\begingroup$ @SolverMax, why not try to solve your problem by LR in an algorithmic manner, e.g. subgradient optimization? If you have still been lucky, you may get the optimal solution. $\endgroup$
    – A.Omidi
    Dec 5, 2022 at 20:14
  • $\begingroup$ I could do that. For example, I have a related model that represents the same situation but assumes Cournot gaming instead of perfect competition. That model can be solved as either an MCP or iteratively as a series of non-linear models that (hopefully) converge to a Nash equilibrium. But my goal here is to solve the perfect competition model using MCP in Pyomo, like I can in GAMS. $\endgroup$
    – Solver Max
    Dec 5, 2022 at 21:17
4
$\begingroup$

You need to get the inequalities correct for each complementarity constraint, and you need to use pyo.inequality to pass a double-sided inequality around the production variables:

import pyomo.environ as pyo
from pyomo import mpec

model = pyo.ConcreteModel()
GENERATORS = [0, 1]
var_cost = [1.0, 1.0]
elasticity, pSet, qSet = -0.3, 20, 3_000
Rho = pSet / elasticity / qSet
model.production = pyo.Var(GENERATORS)
model.demand = pyo.Var()
model.price = pyo.Var()

def c1_rule(model, g):
    return mpec.complements(var_cost[g] - model.price, pyo.inequality(0, model.production[g], 1))

def c2_rule(model):
    return mpec.complements(model.price - (pSet + Rho * (model.demand - qSet)) >= 0, model.demand >= 0)

def c3_rule(model):
    return mpec.complements(sum(model.production[g] for g in GENERATORS) - model.demand == 0, model.price)

model.c1 = mpec.Complementarity(GENERATORS, rule = c1_rule)
model.c2 = mpec.Complementarity(rule = c2_rule)
model.c3 = mpec.Complementarity(rule = c3_rule)

opt = pyo.SolverFactory('pathampl')
results = opt.solve(model, tee = True)

For what it's worth, here's what I get if I solve the problem using JuMP:

using JuMP, PATHSolver
# =====================
# This part is a guess
GENERATORS = 1:2
var_cost = [1.0, 1.0]
g_max = [1.0, 1.0]
# =====================
elasticity, pSet, qSet = -0.3, 20, 3000
ρ = pSet / elasticity / qSet
model = Model(PATHSolver.Optimizer)
@variables(model, begin
    0 <= production[g = GENERATORS] <= g_max[g]
    demand >= 0
    price
end)
@constraints(model, begin
    [g in GENERATORS], var_cost[g] - price ⟂ production[g]
    price - (pSet + ρ * (demand - qSet)) ⟂ demand
    sum(production) - demand ⟂ price
end)
optimize!(model)

julia> solution_summary(model; verbose = true)
* Solver : Path 5.0.03

* Status
  Termination status : LOCALLY_SOLVED
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Result count       : 1
  Has duals          : false
  Message from the solver:
  "The problem was solved"

* Candidate solution
  Objective value      : 0.00000e+00
  Primal solution :
    demand : 2.00000e+00
    price : 8.66222e+01
    production[1] : 1.00000e+00
    production[2] : 1.00000e+00
```
Improve this answer
$\endgroup$
9
  • $\begingroup$ Thanks, with my data your Julia model gets the correct answer. I like Julia, especially how it can succinctly represent a model like this. But changing the Pyomo model to use bounds rather than constraints has no effect - the Pyomo model is still infeasible. $\endgroup$
    – Solver Max
    Dec 6, 2022 at 3:13
  • $\begingroup$ Hello Oscar, any idea if I can run Julia on Google colab, just like python and Gurobi? $\endgroup$
    – Sutanu Majumdar
    Dec 6, 2022 at 3:41
  • 2
    $\begingroup$ I updated my answer with the solution for Pyomo. No idea on Google colab; I haven't used it. But if you google "colab + Julia" you'll find things like this: colab.research.google.com/github/ageron/julia_notebooks/blob/… $\endgroup$
    – Oscar Dowson
    Dec 6, 2022 at 19:30
  • 1
    $\begingroup$ Thanks Oscar. @SolverMax will probably be best to say if it works. I am still novice for both Pyomo & MCP. But it is a pleasure to interact with you. I visited your github and went through Stochastic dual-dynamic programming paper and model. Great contribution to the field. $\endgroup$
    – Sutanu Majumdar
    Dec 6, 2022 at 19:42
  • $\begingroup$ @OscarDowson Using inequality() is a step forward, as it allows the model to solve. However, the solution is wrong. The solution should be: demand = 2100, price = 40 given data: GENERATORS = [0, 1, 2, 3] var_cost = [10.0, 20.0, 30.0, 40.0] GMax = [600.0, 400.0, 1000.0, 2000.0] Both GAMS and Julia find that solution, while the Pyomo solution is: demand = 2637.93, price = 28.05 $\endgroup$
    – Solver Max
    Dec 6, 2022 at 23:27
3
$\begingroup$

In GAMS and Pyomo please check the order of the equations in the complementary constraints lest you are passing $-F$ instead of $F$ where F is the function expressing the constraints. Based on this guide

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks. I've tried several variations, but all produce the same behaviour. $\endgroup$
    – Solver Max
    Dec 5, 2022 at 21:18
1
$\begingroup$

To answer my own question, the following model seems to work as intended.

It seems that GAMS is more forgiving than Pyomo. For the bound on c1_rule's constraint, GAMS allows = or <=, while Pyomo requires <=. In effect, the complementarity relationship says that for a generator's quantity to be non-zero, the generator's cost must be <= the market price. Which is what is intended.

import pyomo.environ as pyo
from pyomo import mpec
import pandas as pd

model = pyo.ConcreteModel()
GENERATORS = [0, 1, 2, 3]
var_cost = [10.0, 20.0, 30.0, 40.0]
GMax = [600.0, 400.0, 1000.0, 2000.0]
elasticity, pSet, qSet = -0.3, 20.0, 3000.0
Rho = pSet / elasticity / qSet

model.production = pyo.Var(GENERATORS, within = pyo.NonNegativeReals, initialize = 0)
model.demand = pyo.Var(within = pyo.NonNegativeReals, initialize = 3000)
model.price = pyo.Var(initialize = 20)

def c1_rule(model, g):
    return mpec.complements(var_cost[g] - model.price <= 0, model.production[g] <= GMax[g])

def c2_rule(model):
    return mpec.complements(model.price - (pSet + Rho * (model.demand - qSet)) == 0, model.demand >= 0)

def c3_rule(model):
    return mpec.complements(sum(model.production[g] for g in GENERATORS) - model.demand == 0, model.price >= 0)

model.c1 = mpec.Complementarity(GENERATORS, rule = c1_rule)
model.c2 = mpec.Complementarity(rule = c2_rule)
model.c3 = mpec.Complementarity(rule = c3_rule)

opt = pyo.SolverFactory('pathampl')
results = opt.solve(model, tee = True)

pd.options.display.float_format = "{:,.2f}".format
GenResults = pd.DataFrame()
for s in GENERATORS:
    GenResults.loc[s, 'Dispatch'] = pyo.value(model.production[s])
    GenResults.loc[s, 'UB'] = GMax[s]
    GenResults.loc[s, 'SRMC'] = var_cost[s]
display(GenResults)

print('Rho:   ', round(Rho,4))
print('Price: ', round(pyo.value(model.price), 2))
print('Demand:', round(pyo.value(model.demand), 2))
Improve this answer
$\endgroup$
2
  • $\begingroup$ I don't think this is the same formulation. For example, I don't know how your third constraint can have both ==0 on the left-hand side and >= 0 on the right-hand side. This is one reason why JuMP does not let people model MCPs with inequalities on both sides. $\endgroup$
    – Oscar Dowson
    Dec 7, 2022 at 4:02
  • $\begingroup$ The effect of having ==0 on the LHS is that the RHS is unconstrained (i.e. it can take any positive or negative value, including zero). But because the production and demand variables are defined as NonNegativeReals, the price cannot be negative. Therefore, in this model, the >= 0 on the RHS doesn't do anything so it can be omitted. The JuMP notation is clearer and more succinct - it looks like how the math formulation would be written, which has a lot of merit. $\endgroup$
    – Solver Max
    Dec 7, 2022 at 7:49
-1
$\begingroup$

In GAMS version I don't see production[g] forming part of the equation unless marginalCost covers it. But you are then setting mpec solver to set either production for g 0 or make marginalCost equal to price.

Improve this answer
$\endgroup$
1
  • $\begingroup$ See my answer, which I think works OK. $\endgroup$
    – Solver Max
    Dec 7, 2022 at 3:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.