python - Pyomo accesing/retrieving dual variables - shadow price with binary variables

Question

I am pretty new to optimization in general and pyomo in particular, so I apologize in advance for any rookie mistakes.

I have defined a simple unit commitment exercise (example 3.1 from [1]) using [2] as starting point. I got the correct result and my code runs, but I have a few questions regarding how to access stuff.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import shutil
import sys
import os.path
import pyomo.environ as pyo
import pyomo.gdp as gdp #necessary if you use booleans to select active and innactive units

def bounds_rule(m, n, param='Cap_MW'):
    # m because it pases the module
    # n because it needs a variable from each set, in this case there was only m.N
    return (0, Gen[n][param]) #returns lower and upper bounds.

def unit_commitment():
    m = pyo.ConcreteModel()
    m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT_EXPORT)
    N=Gen.keys()
    m.N = pyo.Set(initialize=N)    
    
    m.Pgen = pyo.Var(m.N, bounds = bounds_rule) #amount of generation
    m.Rgen = pyo.Var(m.N, bounds = bounds_rule) #amount of generation

    # m.OnOff = pyo.Var(m.N, domain=pyo.Binary) #boolean on/off marker
    
    # objective
    m.cost = pyo.Objective(expr = sum( m.Pgen[n]*Gen[n]['energy_$MWh'] + m.Rgen[n]*Gen[n]['res_$MW'] for n in m.N), sense=pyo.minimize)
    
    # demand
    m.demandP = pyo.Constraint(rule=lambda m: sum(m.Pgen[n] for n in N) == Demand['ener_MWh'])
    m.demandR = pyo.Constraint(rule=lambda m: sum(m.Rgen[n] for n in N) == Demand['res_MW'])
    
    # machine production limits
    # m.lb = pyo.Constraint(m.N, rule=lambda m, n: Gen[n]['Cap_min']*m.OnOff[n] <= m.Pgen[n]+m.Rgen[n] )
    # m.ub = pyo.Constraint(m.N, rule=lambda m, n: Gen[n]['Cap_MW']*m.OnOff[n] >= m.Pgen[n]+m.Rgen[n])
    m.lb = pyo.Constraint(m.N, rule=lambda m, n: Gen[n]['Cap_min'] <= m.Pgen[n]+m.Rgen[n] )
    m.ub = pyo.Constraint(m.N, rule=lambda m, n: Gen[n]['Cap_MW'] >= m.Pgen[n]+m.Rgen[n])

    m.rc = pyo.Suffix(direction=pyo.Suffix.IMPORT)
    
    return m

Gen = {
    'GenA' : {'Cap_MW': 100, 'energy_$MWh': 10, 'res_$MW': 0, 'Cap_min': 0},
    'GenB' : {'Cap_MW': 100, 'energy_$MWh': 30, 'res_$MW': 25, 'Cap_min': 0},
} #starting data

Demand = {
    'ener_MWh': 130, 'res_MW': 20,
} #starting data
m = unit_commitment()
pyo.SolverFactory('glpk').solve(m).write() 
m.display()


df = pd.DataFrame.from_dict([m.Pgen.extract_values(), m.Rgen.extract_values()]).T.rename(columns={0: "P", 1: "R"})
print(df)
print("Cost Function result: " + str(m.cost.expr()) + "$.")

print(m.rc.display())
print(m.dual.display())
print(m.dual[m.demandR])

da= {'duals': m.dual[m.demandP],
     'uslack': m.demandP.uslack(),
     'lslack': m.demandP.lslack(),
    }
db= {'duals': m.dual[m.demandR],
     'uslack': m.demandR.uslack(),
     'lslack': m.demandR.lslack(),
    }

duals = pd.DataFrame.from_dict([da, db]).T.rename(columns={0: "demandP", 1: "demandR"})

print(duals)

Here come my questions.

1-Duals/shadow-price: By definition the shadow price are the dual variables of the constraints (m.demandP and m.demandR). Is there a way to access this values and put them into a dataframe without doing that "shitty" thing I did? I mean defining manually da and db and then creating the dataframe as both dictionaries joined? I would like to do something cleaner like the df that holds the results of P and R for each generator in the system.

2-Usually, in the unit commitment problem, binary variables are used in order to "mark" or "select" active and inactive units. Hence the "m.OnOff" variable (commented line). For what I found in [3], duals don't exist in models containing binary variables. After that I rewrote the problem without including binarys. This is not a problem in this simplistic exercise in which all units run, but for larger ones. I need to be able to let the optimization decide which units will and won't run and I still need the shadow-price. Is there a way to obtain the shadow-price/duals in a problem containing binary variables? I let the constraint definition based on binary variables also there in case someone finds it useful.

Note: The code also runs with the binary variables and gets the correct result, however I couldn't figure out how to get the shadow-price. Hence my question.

[1] Morales, J. M., Conejo, A. J., Madsen, H., Pinson, P., & Zugno, M. (2013). Integrating renewables in electricity markets: operational problems (Vol. 205). Springer Science & Business Media.

[2] https://jckantor.github.io/ND-Pyomo-Cookbook/04.06-Unit-Commitment.html

[3] Dual Variable Returns Nothing in Pyomo

Answer 1

To answer 1, you can dynamically get the constraint objects from your model using model.component_objects(pyo.Constraint) which will return an iterator of your constraints, which keeps your from having to hard-code the constraint names. It gets tricky for indexed variables because you have to do an extra step to get the slacks for each index, not just the constraint object. For the duals, you can iterate over the keys attribute to retrieve those values.

duals_dict = {str(key):m.dual[key] for key in m.dual.keys()}

u_slack_dict = {
    # uslacks for non-indexed constraints
    **{str(con):con.uslack() for con in m.component_objects(pyo.Constraint)
       if not con.is_indexed()},
    # indexed constraint uslack
    # loop through the indexed constraints
    # get all the indices then retrieve the slacks for each index of constraint
    **{k:v for con in m.component_objects(pyo.Constraint) if con.is_indexed()
       for k,v in {'{}[{}]'.format(str(con),key):con[key].uslack()
                   for key in con.keys()}.items()}
    }

l_slack_dict = {
    # lslacks for non-indexed constraints
    **{str(con):con.lslack() for con in m.component_objects(pyo.Constraint)
       if not con.is_indexed()},
    # indexed constraint lslack
    # loop through the indexed constraints
    # get all the indices then retrieve the slacks for each index of constraint
    **{k:v for con in m.component_objects(pyo.Constraint) if con.is_indexed()
       for k,v in {'{}[{}]'.format(str(con),key):con[key].lslack()
                   for key in con.keys()}.items()}
    }

# combine into a single df
df = pd.concat([pd.Series(d,name=name)
           for name,d in {'duals':duals_dict,
                          'uslack':u_slack_dict,
                          'lslack':l_slack_dict}.items()],
          axis='columns')

Regarding 2, I agree with @Erwin s comment about solving with the binary variables to get the optimal solution, then removing the binary restriction but fixing the variables to the optimal values to get some dual variable values.