如何使用 Pyomo 运行这个简单的优化程序？

Question

我正在尝试运行一个“简单”的优化程序，没有任何限制来让程序运行。程序运行后，我计划添加约束和其他输入变量。但是，我似乎无法让该模型发挥作用。为了计算可行性，我已将使用的日期范围过滤到仅 2 天。名为 df 的数据框包含以 15 分钟为间隔的日期（日期时间）以及两列包含建筑负载和光伏发电的数据，均为浮点数（参见附图）。

我使用以下代码只是为了看看优化问题是否有效。请注意，没有添加任何约束（Pgrid_max 被注释掉）。最终我将添加一个包含电动汽车充电数据的数据集，以优化双向电动汽车充电。

start_date = '2023-06-01 00:00:00'
end_date = '2023-06-02 23:59:59'

# Pgrid_max = (df['Building load [kWh]'].max() + 24*18*0.25) * 1.5

model = ConcreteModel()

# Time intervals
T = df[start_date:end_date].index.tolist() # Create a list of the timesteps

# Sets
model.T = Set(initialize=T) # Time interval

# Parameters
model.Pbuilding = Param(model.T, initialize=df[start_date:end_date]['Building load [kWh]'].to_dict())
model.Ppv = Param(model.T, initialize=df[start_date:end_date]['PV load [kWh]'].to_dict())


# Variables
model.Pgrid = Var(model.T, domain=Reals)
# model.Pev = Var(model.T, domain=PositiveReals, bounds=(0, 24*0.25)) 


# Energy balance equation
def energy_balance(model, t):
    return model.Pgrid[t] == model.Pbuilding[t] - model.Ppv[t]

# Objective function
def objective_rule(model):
    return sum(model.Pgrid[t] for t in model.T)

model.obj = Objective(rule=objective_rule, sense=minimize)

# Solve the model (assuming you have a solver like glpk installed)
solver = SolverFactory('glpk')
solver.solve(model, tee=True)

# Extract results directly into a pandas DataFrame
results = pd.DataFrame(index=T)
results['Pgrid'] = [model.Pgrid[t].value for t in T]

for t in model.T:
    results['Pgrid'][t] = model.Pgrid[t].value

print(results)

运行此代码时，我不断收到问题无法解决的消息。 Pgrid 可以取任何实际值吗？打印结果仅给出“无”值。有人知道我做错了什么吗？

Answer 1

您的代码中有几个问题。

最重要的（也是为什么你可能得到零）是你没有构建约束。您需要将该函数作为

rule

调用，如下面的代码所示。

其他一些技巧/想法：

• 始终

  pprint

  对模型进行质量检查。如果你这样做了，你会注意到缺少的约束
• 始终检查以确保结果是最佳的，否则结果就是垃圾。您可以检查求解器的措辞或放入断言，如下所示
• 当事情变得混乱时，请将

  pandas

  收起来，除非您对此感到非常舒服。如图所示的字典很容易用小数据排除故障
• 如果您打算在完成后将内容放回 DF 中，我建议您像我展示您在哪里使用字典一样进行操作，以确保索引正确。

在您的型号上：

• 您可能希望将

  Pgrid

  设置为非负实数，并使余额约束大于或等于 (GTE)，如图所示，以涵盖 pv > 需求的情况，除非您可以将功率推回至网格。

代码：

"""
Simple Power Grid Model
"""
from pyomo.environ import *
import pandas as pd

### DATA

# take pandas out of the equation for simple model...
load = { 0:  8.2,
         15: 9.5,
         30: 7.7,
         45: 4.3}

pv = {   0:  4.2,
         15: 5.1,
         30: 6.0,
         45: 8.8}

# start_date = '2023-06-01 00:00:00'
# end_date = '2023-06-02 23:59:59'

model = ConcreteModel()

# Time intervals
# T = df[start_date:end_date].index.tolist() # Create a list of the timesteps

# Sets
model.T = Set(initialize=load.keys()) # Time interval

# Parameters
model.Pbuilding = Param(model.T, initialize=load) # df[start_date:end_date]['Building load [kWh]'].to_dict())
model.Ppv = Param(model.T, initialize=pv) # df[start_date:end_date]['PV load [kWh]'].to_dict())


# Variables
model.Pgrid = Var(model.T, domain=NonNegativeReals) # probably can't "push" back to grid...
# model.Pev = Var(model.T, domain=PositiveReals, bounds=(0, 24*0.25))


# Energy balance equation
def energy_balance(model, t):
    # make this a GTE constraint:  must pull enough from grid to make up deficit, if there is one
    return model.Pgrid[t] >= model.Pbuilding[t] - model.Ppv[t]
# you were MISSING the contraint construction call w/ a rule:
model.balance = Constraint(model.T, rule=energy_balance)

# Objective function
def objective_rule(model):
    return sum(model.Pgrid[t] for t in model.T)
model.obj = Objective(rule=objective_rule, sense=minimize)

#  INSPECT THE MODEL!!!
model.pprint()

# Solve the model (assuming you have a solver like glpk installed)
solver = SolverFactory('cbc')  # or glpk or ...
results = solver.solve(model, tee=True)  # catch the results object

#  Ensure that it was OPTIMAL or results are junk
status = results.Solver()['Termination condition'].value
assert status == 'optimal', f'error occurred, status: {status}.  Check model!'

# Extract results directly into a pandas DataFrame
results = pd.DataFrame(index=model.T)

# providing the item pairs is "safer" as there is no risk of issues
# with the order of the values pulled from the set model.T
results['Pgrid'] = pd.Series({k:v.value for k, v in model.Pgrid.items()}) # [model.Pgrid[t].value for t in model.T]

print(results)

输出：

1 Set Declarations
    T : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    4 : {0, 15, 30, 45}

2 Param Declarations
    Pbuilding : Size=4, Index=T, Domain=Any, Default=None, Mutable=False
        Key : Value
          0 :   8.2
         15 :   9.5
         30 :   7.7
         45 :   4.3
    Ppv : Size=4, Index=T, Domain=Any, Default=None, Mutable=False
        Key : Value
          0 :   4.2
         15 :   5.1
         30 :   6.0
         45 :   8.8

1 Var Declarations
    Pgrid : Size=4, Index=T
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  None :  None : False :  True : NonNegativeReals
         15 :     0 :  None :  None : False :  True : NonNegativeReals
         30 :     0 :  None :  None : False :  True : NonNegativeReals
         45 :     0 :  None :  None : False :  True : NonNegativeReals

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : Pgrid[0] + Pgrid[15] + Pgrid[30] + Pgrid[45]

1 Constraint Declarations
    balance : Size=4, Index=T, Active=True
        Key : Lower              : Body      : Upper : Active
          0 :  3.999999999999999 :  Pgrid[0] :  +Inf :   True
         15 :                4.4 : Pgrid[15] :  +Inf :   True
         30 : 1.7000000000000002 : Pgrid[30] :  +Inf :   True
         45 : -4.500000000000001 : Pgrid[45] :  +Inf :   True

6 Declarations: T Pbuilding Ppv Pgrid balance obj
Welcome to the CBC MILP Solver 
Version: 2.10.10 
Build Date: Jul 24 2023 

command line - /opt/homebrew/opt/cbc/bin/cbc -printingOptions all -import /var/folders/7l/f196n6c974x3yjx5s37t69dc0000gn/T/tmp28wkyeg0.pyomo.lp -stat=1 -solve -solu /var/folders/7l/f196n6c974x3yjx5s37t69dc0000gn/T/tmp28wkyeg0.pyomo.soln (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 0 (-4) rows, 0 (-4) columns and 0 (-4) elements
Statistics for presolved model


Problem has 0 rows, 0 columns (0 with objective) and 0 elements
Column breakdown:
0 of type 0.0->inf, 0 of type 0.0->up, 0 of type lo->inf, 
0 of type lo->up, 0 of type free, 0 of type fixed, 
0 of type -inf->0.0, 0 of type -inf->up, 0 of type 0.0->1.0 
Row breakdown:
0 of type E 0.0, 0 of type E 1.0, 0 of type E -1.0, 
0 of type E other, 0 of type G 0.0, 0 of type G 1.0, 
0 of type G other, 0 of type L 0.0, 0 of type L 1.0, 
0 of type L other, 0 of type Range 0.0->1.0, 0 of type Range other, 
0 of type Free 
Presolve 0 (-4) rows, 0 (-4) columns and 0 (-4) elements
Empty problem - 0 rows, 0 columns and 0 elements
Optimal - objective value 10.1
After Postsolve, objective 10.1, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 10.1 - 0 iterations time 0.002, Presolve 0.00
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

    Pgrid
0     4.0
15    4.4
30    1.7
45    0.0

Process finished with exit code 0