如何在PYOMO的目标函数中同时添加参数值和变量值？

Question

我写了一个程序，目的是为了最大化收入，这个程序涉及一个设施，里面有太阳能和风能发电厂。它们的任务是满足一定的需求。为了满足这个需求，它们会根据PPA（电力购买协议）获得一些收入。影响收入的因素有三个：

1. 实际收入：最小值（设施发电量（太阳能 + 风能），需求） * PPA
2. 额外收入：（最大值（设施发电量，需求） - 需求） * PPA的50%
3. 短缺罚款（如果设施的输出低于需求的90%（DFR），那么就会适用这个）：（最小值（设施输出，DFR * 需求） - 需求 * 短缺罚款

总收入 = 实际收入 + 额外收入 - 短缺罚款 这是我为此写的程序：

import numpy as np
import pandas as pd
from pyomo.environ import *

total_instances = 8760
np.random.seed(total_instances)
load_profile = np.random.randint(80, 130, total_instances)
solar_profile = np.random.uniform(0, 0.9, total_instances)
wind_profile = np.random.uniform(0, 0.9, total_instances)
month_idx = pd.date_range('1/1/2023', periods=total_instances, freq='60min').month.tolist()
PPA = 10
shortfall_penalty_rate = 15
excess_rate = 5
DFR = 0.9
solar_capacity = 120
wind_capacity = 130

# Not in Use right now
load_profile_month = list(zip(month_idx, load_profile))
monthly_sum = [0] * 12
for i in range(len(load_profile_month)):
    month_idx, load_value = load_profile_month[i]
    monthly_sum[month_idx - 1] += load_value

model = ConcreteModel()

model.m_index = Set(initialize=list(range(len(load_profile))))

# variable
model.facility_output = Var(model.m_index, domain=NonNegativeReals)

#  gen variable
model.solar_use = Var(model.m_index, domain=NonNegativeReals)
model.wind_use = Var(model.m_index, domain=NonNegativeReals)

# Load profile
model.load_profile = Param(model.m_index, initialize=load_profile)
model.solar_profile = Param(model.m_index, initialize=solar_profile)
model.wind_profile = Param(model.m_index, initialize=wind_profile)

model.lost_load = Var(model.m_index, domain=NonNegativeReals)


# Objective function
def revenue(model):
    actual_revenue = sum(
        min((model.solar_use[m] + model.wind_use[m]), model.load_profile[m]) * PPA
        for m in model.m_index
    )
    excess_revenue = sum(
        (max(model.solar_use[m] + model.wind_use[m], model.load_profile[m]) - model.load_profile[m]) * excess_rate
        for m in model.m_index
    )
    shortfall_penalty = sum(
        (min(model.solar_use[m] + model.wind_use[m], DFR * model.load_profile[m]) - model.load_profile[m] * DFR) * shortfall_penalty_rate
        for m in model.m_index
    )
    total_revenue = actual_revenue + excess_revenue + shortfall_penalty

    return total_revenue


model.obj = Objective(rule=revenue, sense=maximize)


def energy_balance(model, m):
    return model.grid[m] == model.solar_use[m] + model.wind_use[m] + model.lost_load[m]


model.energy_balance = Constraint(model.m_index, rule=energy_balance)


def grid_limit(model, m):
    return model.grid[m] >= model.load_profile[m]


model.grid_limit = Constraint(model.m_index, rule=grid_limit)


def max_solar(model, m):
    eq = model.solar_use[m] <= solar_capacity * model.solar_profile[m]
    return eq


model.max_solar = Constraint(model.m_index, rule=max_solar)


def max_wind(model, m):
    eq = model.wind_use[m] <= wind_capacity * model.wind_profile[m]
    return eq


model.max_wind = Constraint(model.m_index, rule=max_wind)

Solver = SolverFactory('gurobi')
Solver.options['LogFile'] = "gurobiLog"
# Solver.options['MIPGap'] = 0.0
print('\nConnecting to Gurobi Server...')
results = Solver.solve(model)

if (results.solver.status == SolverStatus.ok):
    if (results.solver.termination_condition == TerminationCondition.optimal):
        print("\n\n***Optimal solution found***")
        print('obj returned:', round(value(model.obj), 2))
    else:
        print("\n\n***No optimal solution found***")
        if (results.solver.termination_condition == TerminationCondition.infeasible):
            print("Infeasible solution")
            exit()
else:
    print("\n\n***Solver terminated abnormally***")
    exit()

grid_use = []
solar = []
wind = []
# e_in = []
# e_out = []
# soc = []
lost_load = []
load = []
for i in range(len(load_profile)):
    grid_use.append(value(model.grid[i]))
    solar.append(value(model.solar_use[i]))
    wind.append(value(model.wind_use[i]))
    lost_load.append(value(model.lost_load[i]))
    load.append(value(model.load_profile[i]))
    # e_in.append(value(model.e_in[i]))
    # e_out.append(value(model.e_out[i]))
    # soc.append(value(model.soc[i]))
pd.DataFrame({
    'Grid': grid_use,
    'Solar': solar,
    'Wind': wind,
    'Shortfall': lost_load,
    'Load Profile': load,
}).to_excel('testing4.xlsx')

我遇到的错误，我认为是因为我在使用负载配置，它是一个参数，与变量中的太阳能和风能有关：

ValueError: Error retrieving immutable Param value (load_profile[0]):
    The Param value is undefined and no default value is specified.
ERROR: Rule failed when generating expression for objective obj: ValueError:
    Error retrieving immutable Param value (load_profile[0]):
        The Param value is undefined and no default value is specified.
ERROR: Constructing component 'obj' from data=None failed: ValueError: Error
    retrieving immutable Param value (load_profile[0]):
        The Param value is undefined and no default value is specified.

Process finished with exit code 1

Answer 1

你用来初始化参数的所有数据项应该是字典（或者其他数据类型，比如pandas的序列），这些字典里要有键值对。你的参数是有索引的，所以你输入的数据需要用相同的键来标识。

favorability = {'dog': 10, 'cat': 2, 'pig': 4}
random_numbers = { 1: 22.5, 2: 52.9, 3: 4.1}

在模型中，类似这样的：

m.animals = m.Set(initialize=favorability.keys())
m.case = m.Set(initialize=random_numbers.keys())

m.favorability = m.Param(m.animals, initialize=favorability)
m.cost = m.Param(m.case, initialize=random_numbers)