Search code examples
pythonoptimizationpyomo

Pyomo Building a Lagrangian or Selecting left hand side of constraints

I have a pyomo model "m" with 4 variables and several constraints (both equality and inequality) in the form:

Min F(G1,G2,D1,D2) st h=0 g<=0

Then I need to build the lagrangian function, which is something like this: enter image description here

Briefly, lambda and mu are the duals of the constraints. So I need the objective function + dual1cons1 + dual2cons2 and so on.

I have literally no idea how to do this. The closest I got is with this:

Lagrange = m.objective #initialize the variable with the f(x)

for cons in m.component_objects(pyomo.core.base.constraint.Constraint): #iterates over the constraints
    print("\n->" + str(cons) + ":")
    if isinstance(cons, pyomo.core.base.constraint.SimpleConstraint): #selects whether it is an individual constraint
        # print(cons)
        print(expression_to_string(m.component(cons).expr)) #prints the expresion
    if isinstance(cons, pyomo.core.base.constraint.ConstraintList): #or a list of constraints
        for i in m.component(cons):
            Lagrange=Lagrange+m.component('LinLim')[i].expr
            # print(expression_to_string(m.component(cons)[i].expr)) #prints the expresion

print(expression_to_string(Lagrange))

Then, the print statement returns this:

(12Pgen[G1] + 20Pgen[G2] - (40Pdem[D1] + 35Pdem[D2])) + (500*(teta[n1] - teta[n2]) - 100 <= 0.0) - (500*(teta[n1] - teta[n2]) + 100) <= 0.0 + (500*(teta[n1] - teta[n3]) - 100 <= 0.0) + (500*(teta[n1] - teta[n2]) - 100 <= 0.0) - (500*(teta[n1] - teta[n2]) + 100) <= 0.0 + (500*(teta[n1] - teta[n3]) - 100 <= 0.0) - (500*(teta[n1] - teta[n3]) + 100) <= 0.0 + (500*(teta[n2] - teta[n3]) - 100 <= 0.0) - (500*(teta[n2] - teta[n3]) + 100) <= 0.0

I am aware that is a nightmare to read. The point is, it includes the operator of the equations (==, >=, <=) while I would only be interested in the left hand side part of the equation. Then, in addition I would need to add new variables representing the duals (lambda and mu).

Solution

  • Pyomo constraints have a body attribute, and uslack and lslack functions that give you a sum expression of the constraint left hand side, the upper slack, and lower slack respectively. The `body attribute is what you want to multiply by lambda and mu. Here is an example with a simple constraint

    import pyomo.environ as pyo

m = pyo.ConcreteModel()

index = [1,2,3,4]
m.x = pyo.Var(index,domain = pyo.Binary,initialize = 0)

m.c1 = pyo.Constraint(expr=sum([m.x[i] for i in index])<=3)
print('Constraint 1: {}'.format(str(m.c1.body)))
# Constraint 1: x[1] + x[2] + x[3] + x[4]
m.fixed_lambda = pyo.Var(domain = pyo.NonNegativeReals)
m.lagrangian = pyo.Objective(expr = m.fixed_lambda * m.c1.body,sense=pyo.minimize)
print('Lagrangian expression: {}'.format(str(m.lagrangian.expr)))
# Lagrangian expression: fixed_lambda*(x[1] + x[2] + x[3] + x[4])