Pyomo 访问/检索双变量 - 带二元变量的影子价格 [英] Pyomo accesing/retrieving dual variables - shadow price with binary variables

问题描述

我对一般优化特别是 pyomo 还很陌生，因此对于任何新手错误，我提前道歉.

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

我使用 [2] 作为起点定义了一个简单的单元承诺练习(来自 [1] 的示例 3.1).我得到了正确的结果并且我的代码运行了，但我有一些关于如何访问内容的问题.

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)

我的问题来了.

1-Duals/shadow-price:根据定义，影子价格是约束的对偶变量(m.demandP 和 m.demandR).有没有一种方法可以访问这些值并将它们放入数据帧中而不用做糟糕"的事情?我做了什么?我的意思是手动定义 da 和 db，然后在两个字典加入时创建数据框?我想做一些更清晰的事情，比如保存系统中每个生成器的 P 和 R 结果的 df.

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-通常，在单元提交问题中，使用二元变量来标记"或选择"活动和非活动单位.因此，m.OnOff"变量(注释行).对于我在 [3] 中发现的内容，包含二元变量的模型中不存在对偶.之后，我重新编写了问题而不包含二进制文件.在这个所有单元都运行的简单练习中，这不是问题，但对于更大的单元而言，这不是问题.我需要能够让优化决定哪些单位将运行和不运行，我仍然需要影子价格.有没有办法在包含二元变量的问题中获得影子价格/对偶?我让基于二元变量的约束定义也在那里，以防有人发现它有用.

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).将可再生能源纳入电力市场:运营问题(第 205 卷).施普林格科学与商业媒体.

[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] 双变量在 Pyomo 中没有任何回报

推荐答案

要回答 1，您可以使用 model.component_objects(pyo.Constraint) 从模型中动态获取约束对象，这将返回约束的迭代器，它使您不必对约束名称进行硬编码.索引变量变得棘手，因为您必须执行额外的步骤来获取每个索引的松弛量，而不仅仅是约束对象.对于对偶，您可以遍历 keys 属性以检索这些值.

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')

关于 2，我同意 @Erwin 关于用二进制变量求解以获得最优解，然后去除二进制限制但将变量固定为最优值以获得一些对偶变量值的评论.

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.

这篇关于Pyomo 访问/检索双变量 - 带二元变量的影子价格的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！