如何应用约束以使生成器运行至少 4 个间隔？
问题描述 投票：0回答：1
我正在使用 Pyomo 编写优化问题陈述。这是一个一般性的问题陈述：
- 需求概况为96 个区块。
- 有2台700kW主发电机（1代和2代）和1台1000kW备用发电机。
- 发电机有两个单独的价目表价格 1 和价格 2，备用发电机是最昂贵的。
- 发电机本质上是排他性的，这意味着一次只有一台发电机处于活动状态，其余的能量将来自备用发电机。
- 最后也是最重要的事情是，如果发电机处于活动状态，那么它应该在接下来的 4 个区块中保持活动状态。这就是我面临的问题。
我在 Pyomo 中为此创建了一个具体模型。这是它的代码。
import numpy as np
from pyomo.environ import *
from pyomo.opt import SolverFactory
import pandas as pd
np.random.seed(96)
prices = {
'price 1': np.random.randint(1,20,96),
'price 2': np.random.randint(1,20,96),
}
load_profile = np.random.randint(500,1000,96)
generator_capacities = {
'gen 1': 700,
'gen 2': 700,
'backup gen': 1000,
}
backup_cost = 10000
model = ConcreteModel()
# Define decision variables
model.m_index = Set(initialize=list(range(len(load_profile))))
model.gen1_output = Var(model.m_index, domain=NonNegativeReals)
model.gen2_output = Var(model.m_index, domain=NonNegativeReals)
model.backup_gen_output = Var(model.m_index, domain=NonNegativeReals)
# Define binary variables for gen1 and gen2
model.gen1_active = Var(model.m_index, within=Binary, initialize=1)
model.gen2_active = Var(model.m_index, within=Binary, initialize=0)
model.gen1_continuous_duration = Var(model.m_index, within=Binary)
model.gen2_continuous_duration = Var(model.m_index, within=Binary)
model.load_data = Param(model.m_index, initialize=dict(zip(model.m_index, load_profile)))
def cost(model):
gen1_cost = sum(model.gen1_output[m] * prices['price 1'][m] for m in model.m_index)
gen2_cost = sum(model.gen2_output[m] * prices['price 2'][m] for m in model.m_index)
lost_load = sum(model.backup_gen_output[m] for m in model.m_index)
total_cost = gen1_cost + gen2_cost + lost_load * backup_cost
return total_cost
model.obj = Objective(rule=cost, sense=minimize)
def load_constraint_rule(model, m):
return (
model.gen1_output[m] + model.gen2_output[m] + model.backup_gen_output[m] == load_profile[m]
)
model.load_constraint = Constraint(model.m_index, rule=load_constraint_rule)
def gen1_max(model, m):
eq = model.gen1_output[m] <= 700 * model.gen1_active[m]
return eq
model.gen1_max = Constraint(model.m_index, rule=gen1_max)
def gen2_max(model, m):
eq = model.gen2_output[m] <= 700 * model.gen2_active[m]
return eq
model.gen2_max = Constraint(model.m_index, rule=gen2_max)
def backup_max(model, m):
eq = model.backup_gen_output[m] <= 1000
return eq
model.backup_max = Constraint(model.m_index, rule=backup_max)
def exclusive_gen_operation_rule(model, m):
return model.gen1_active[m] + model.gen2_active[m] <= 1
model.exclusive_gen_operation = Constraint(model.m_index, rule=exclusive_gen_operation_rule)
def gen1_continuous(model, m):
interval = 4
if m == 0:
return Constraint.Skip
elif m == 1:
if (model.gen1_active[m-1] == 1):
if m <= 96 - interval:
eq = sum([model.gen1_active[i] for i in range(m, m + interval)]) == interval
return eq
else:
return Constraint.Skip
else:
return Constraint.Skip
elif m >= 2:
if (model.gen1_active[m-1] == 1 and model.gen1_active[m-2] == 0):
if m <= 96 - interval:
eq = sum([model.gen1_active[i] for i in range(m, m + interval)]) == interval
return eq
else:
return Constraint.Skip
else:
return Constraint.Skip
model.gen1_continuous = Constraint(model.m_index, rule=gen1_continuous)
def gen2_continuous(model, m):
interval = 4
if m == 0:
return Constraint.Skip
elif m == 1:
if (model.gen2_active[m-1] == 1):
if m <= 96 - interval:
eq = sum([model.gen2_active[i] for i in range(m, m + interval)]) == interval
return eq
else:
return Constraint.Skip
else:
return Constraint.Skip
elif m >= 2:
if (model.gen2_active[m-1] == 1 and model.gen2_active[m-2] == 0):
if m <= 96 - interval:
eq = sum([model.gen2_active[i] for i in range(m, m + interval)]) == interval
return eq
else:
return Constraint.Skip
else:
return Constraint.Skip
model.gen2_continuous = Constraint(model.m_index, rule=gen2_continuous)
# def duration_constraint1_rule(model, m):
# T = 4 # Desired minimum duration
# if m <= 96 - T:
# return sum(model.gen1_active[i] for i in range(m, m + T)) == T
# return Constraint.Skip
#
# model.duration_constraint1 = Constraint(model.m_index, rule=duration_constraint1_rule)
#
# def duration_constraint2_rule(model, m):
# T = 4 # Desired minimum duration
# if m <= 96 - T:
# return sum(model.gen2_active[i] for i in range(m, m + T)) == T
# return Constraint.Skip
#
# model.duration_constraint2 = Constraint(model.m_index, rule=duration_constraint2_rule)
#
#
# def continuous_sync_rule(model, m):
# return model.gen1_active[m] == model.gen1_continuous_duration[m]
#
# model.continuous_sync1 = Constraint(model.m_index, rule=continuous_sync_rule)
#
# def continuous_sync_rule(model, m):
# return model.gen2_active[m] == model.gen2_continuous_duration[m]
#
# model.continuous_sync2 = Constraint(model.m_index, rule=continuous_sync_rule)
opt = SolverFactory('gurobi') # Use a solver of your choice
results = opt.solve(model)
if results.solver.termination_condition == TerminationCondition.optimal:
print("Optimal solution found")
else:
print("Solver did not find an optimal solution")
output={}
for gen in ['gen1', 'gen2', 'backup_gen']:
values = [model.__getattribute__(gen + '_output')[hour].value for hour in model.m_index]
output[gen] = values
# print(f"{gen} output: {values}")
for gen in ['gen1', 'gen2']:
values = [model.__getattribute__(gen + '_active')[hour].value for hour in model.m_index]
output[gen + ' switch'] = values
# print(f"{gen} active:", values)
output['price 1'] = prices['price 1']
output['price 2'] = prices['price 2']
output['Load profile'] = load_profile
pd.DataFrame(output).to_excel('gen switches.xlsx')
# print(prices)
gen1_连续和gen2_连续是我尝试应用此条件的约束，如果gen1处于活动状态，那么它应该在接下来的4个块中保持活动状态，但它无法正常工作。有人可以帮忙吗？
预期的输出应该是这样的。这是模型的输出（因为我使用随机数据，你可能会得到不同的输出）：
但是，根据您在上图中看到的价格，gen1 switch和gen2 switch的预期输出应该是这样的。
1个回答
0
投票
投票
这是我认为有效并能解决您的限制的策略。正如评论中提到的，我认为使用两个具有链接的变量来建模“启动”和“运行”之间的关系是有效的
Start => forces next 4 (or till end) periods to be running
Running <= only if started or running in previous period
双重索引这也是谨慎的做法，特别是当您要添加更多源时。我用你的数据移动了一些小家具，并稍微降低了最小负载，因为备份似乎总是在原始数据中运行，因为在任何 4 个周期中负载高于 700 且至少为 500 的可能性相当高。
我在目标中加了一点“糖”以停止不必要的运行。如果您删除它，备份将在第一个周期开始，并在整个时间内愉快地闲置，在必要时增加负载。可能不想要。
代码：
# Generation Scheduling
import numpy as np
from pyomo.core import sum_product
from pyomo.environ import Var, Set, NonNegativeReals, Binary, Param, Objective, minimize, ConcreteModel, Constraint, \
SolverFactory, TerminationCondition
# from pyomo.opt import SolverFactory
import pandas as pd
# np.random.seed(96)
time_periods = 100
backup_cost = 10000
prices = {
(gen, time): np.random.randint(1, 20) for gen in ['Bertha', 'Sammy'] for time in range(time_periods)
}
for t in range(time_periods):
prices['Backup', t] = backup_cost
load_profile = np.random.randint(200, 1000, time_periods)
generator_capacities = {
'Bertha': 700,
'Sammy': 700,
'Backup': 1000,
}
model = ConcreteModel()
### SETS
model.T = Set(initialize=list(range(len(load_profile))), doc='time')
model.G = Set(initialize=generator_capacities.keys(), doc='generator')
### VARS
model.output = Var(model.G, model.T, domain=NonNegativeReals)
# Define binary variables for gen1 and gen2
model.gen_active = Var(model.G, model.T, domain=Binary)
model.gen_start = Var(model.G, model.T, domain=Binary)
### PARAMS
model.demand = Param(model.T, initialize=dict(zip(model.T, load_profile)))
model.capacity = Param(model.G, initialize=generator_capacities)
model.price = Param(model.G, model.T, initialize=prices)
### OBJ
# cost to idle generator. If you don't do this, the backup will always be "active" even if no output
active_penalty = 0.1
model.obj = Objective(expr=sum_product(model.price, model.output) +
active_penalty * sum(model.gen_active[g, t] for g in model.G for t in model.T))
### Constraints
@model.Constraint(model.T)
def meet_demand(model, t):
return sum(model.output[g, t] for g in model.G) >= model.demand[t]
@model.Constraint(model.G, model.T)
def dont_overspeed_generator(model, g, t):
return model.output[g, t] <= model.capacity[g] * model.gen_active[g, t]
@model.Constraint(model.T)
def select_generator(model, t):
return sum(model.gen_active[g, t] for g in model.G - {'Backup',}) <= 1
@model.Constraint(model.G, model.T)
def min_runtime(model, g, t):
"""any generator that is started must run for 4 periods, unless in final 3 periods"""
run_times = set(range(t, min(t+4, model.T.last())))
return sum(model.gen_active[g, tt] for tt in run_times) >= len(run_times) * model.gen_start[g, t]
@model.Constraint(model.G, model.T)
def continuity(model, g, t):
"""can only be running if started in this period or running in previous"""
if t == model.T.first():
return model.gen_active[g, t] <= model.gen_start[g, t]
else:
return model.gen_active[g, t] <= model.gen_start[g, t] + model.gen_active[g, t-1]
# model.pprint()
opt = SolverFactory('cbc') # Use a solver of your choice
results = opt.solve(model)
if results.solver.termination_condition == TerminationCondition.optimal:
print("Optimal solution found")
else:
print("Solver did not find an optimal solution")
output = {}
output['Load profile'] = load_profile
for gen in model.G:
values = [model.output[gen, hour].value for hour in model.T]
output[f'{gen} output'] = values
# print(f"{gen} output: {values}")
for gen in model.G:
values = [model.gen_active[gen, hour].value for hour in model.T]
output[f'{gen} active'] = values
# print(f"{gen} active:", values)
for gen in model.G :
values = [model.gen_start[gen, hour].value for hour in model.T]
output[f'{gen} start'] = values
for gen in model.G - {'Backup'}:
price_list = []
for t in range(time_periods):
price_list.append(prices[gen, t])
output[f'{gen} price'] = price_list
df = pd.DataFrame(output)
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', None)
pd.set_option('display.width', 300)
print(df)
输出（非常宽）
Load profile Bertha output Sammy output Backup output Bertha active Sammy active Backup active Bertha start Sammy start Backup start Bertha price Sammy price
0 405 405.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 2 11
1 780 700.0 0.0 80.0 1.0 0.0 1.0 0.0 0.0 1.0 19 17
2 620 620.0 0.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0 6 12
3 273 273.0 0.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0 17 7
4 958 700.0 0.0 258.0 1.0 0.0 1.0 0.0 0.0 0.0 8 17
5 777 700.0 0.0 77.0 1.0 0.0 1.0 0.0 0.0 0.0 9 16
6 510 0.0 510.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 14 1
最新问题
- 不确定为什么脚本运行但不给出输出？
- 手动计算TPR、FPR与通过ROC曲线的scikit-learn计算
- 如何修复 Fivem 中的“脚本错误：资源 es_extended 中没有此类导出 getSharedObject”错误？
- onChange(e) 拒绝触发
- 如何在 Gio.Settings 中获取可重定位架构的路径？
- 将每个包含阿拉伯字母的单词更改为反转
- Monorepo-Yarn 工作区“找不到模块或其相应的类型声明。”构建包后
- “JwtBearerDefaults.AuthenticationScheme”和“JwtBearerDefaults”到底是什么？
- 如何强制特定值出现在 SQL 排序中的顶部？
- 找到 ( ) 之间的所有内容并为其创建超链接
- PHP strpos() 函数返回错误结果[重复]
- strpos() 意外返回 false
- 可重用的 React hook 与复杂的 javascript 响应结构
- 从 R 中的键：值对重命名变量
- 如何通过按转义键来转义读取主机？
- 如何处理损坏的 Git 对象文件？
- Zsh | `酿造清单| {任何命令}`引发未捕获的信号管道错误
- 如何对两个列表的元素进行随机分组，并使用随机获得的新元素的元组创建新列表？
- 将 MySQL 数据库导入 SQL Server
- 在 SQL Server 中生成中文的 UTF-8 MD5 哈希值
© www.soinside.com 2019 - 2024. All rights reserved.