Pyomo 优化不起作用（天然气厂调度）
Pyomo optimisation not working (gas plant dispatch)
背景
我正在尝试编写一个 pyomo 脚本，以根据对电价的完美预见来优化调度天然气厂。我相信我已经完成了 90%，只有几个问题。
问题
我的脚本有效，但求解器从不调度工厂，即使它应该在的地方，在下面提供的示例中，我可以手动计算至少 8131 美元的潜在利润。
我怀疑结果为零的原因是我编写约束的方式造成的，其中有 2 个；
- Gas Plant 从冷启动开始需要 10 分钟
- 预热后，燃气设备有一个必须运行的最小负载 at/above。
- 天然气厂一天只能消耗 9000 GJ 的天然气
特别是在进一步测试中，我认为是 'gas_volume_used' 约束导致了问题。
请求帮助
有人可以看看我的代码，看看我在约束方程中遗漏了什么吗？
代码
import pandas as pd
from pyomo.environ import *
from pyomo.core import *
# ---------------------
# Main Functions
# ---------------------
def model_to_df(model, first_period, last_period):
# Need to increase the first & last period by 1 because of pyomo indexing
# and add 1 to the value of last model period because of the range
periods = range(model.T[first_period + 1], model.T[last_period + 1] + 1)
rrp = [model.P.extract_values()[None][i] for i in periods]
Gascold = [value(model.Gascold[i]) for i in periods]
Gaswarm = [value(model.Gaswarm[i]) for i in periods]
Gashot = [value(model.Gashot[i]) for i in periods]
GasUnit_mw = [value(model.GasUnit_mw[i]) for i in periods]
GasUse = [value(model.GasUse[i]) for i in periods]
df_dict = dict(
period=periods,
rrp=rrp,
Gascold=Gascold,
Gaswarm=Gaswarm,
Gashot=Gashot,
GasUnit_mw=GasUnit_mw,
GasUse=GasUse
)
df = pd.DataFrame(df_dict)
return df
def optimize_year(df, gunits, gnet, gprice, vom, cold_gas_p_unit, hr, max_gas, sys_use,
first_model_period, last_model_period):
"""
Optimize behavior of a gas plant over a time period
Assume perfect foresight of electricity prices.
Parameters
----------
df : dataframe
dataframe with columns of hourly LBMP and the hour of the year
gunits : int
number of gas turbines
gnet : double
size in MW's of each gas turbine
gprice : double
the price of gas $/GJ
vom : double
the variable cost of gas $/MWh dispatched
cold_gas_p_unit : double
the amount of gas in GJ used in 5 minutes while booting up from cold start
hr : double
heat rate, the amount of gas consumed per MWh produced (GJ/MWh)
sys_use : double
auxillary gas usage, percentage
max_gas : double
the maximum amount of gas (GJ) available to use in a period
first_model_period : int, optional
Set the first period of the year to be considered in the optimization
last_model_period : int, optional
Set the last period of the year to be considered in the optimization
Returns
-------
dataframe
Gas plant operational stats and time stamp
"""
# Filter the data
df = df.loc[first_model_period:last_model_period, :]
model = ConcreteModel()
# Define model parameters
model.T = Set(initialize=df.period.tolist(), doc='periods of year', ordered=True)
model.P = Param(initialize=df.rrp.tolist(), doc='price for each period', within=Any)
model.MaxGas = Param(initialize=(gnet/12) * gunits * hr * (1+sys_use), doc='Max gas usage (GJ) in interval', within=Any)
model.MaxVol = Param(initialize=max_gas, doc='Max gas usage (GJ) in 24 hour period', within=Any)
model.MaxMW = Param(initialize=(gnet/12) * gunits, doc='Max MWs', within=Any)
model.Gascold = Var(model.T, doc='binary', within=Binary, initialize=1)
model.Gaswarm = Var(model.T, doc='binary', within=Binary, initialize=0)
model.Gashot = Var(model.T, doc='binary', within=Binary, initialize=0)
model.GasUse = Var(model.T, bounds=(0, model.MaxGas))
model.GasUnit_mw = Var(model.T, bounds=(0, model.MaxMW))
# *********************
# GAS CONSTRAINTS
# *********************
"Ensure that all three Gas State flags sum to 1"
def set_on(model, t):
return model.Gascold[t] + model.Gaswarm[t] + model.Gashot[t] == 1
model.limit_flags = Constraint(model.T, rule=set_on)
"Set t0 to GAS COLD = 1"
def gas_cold(model, t):
if t == model.T.first():
return model.Gascold[t] == 1
else:
return model.Gascold[t] >= 0
model.gas_cold_check = Constraint(model.T, rule=gas_cold)
"Gas can only warm up from cold"
def gas_warm(model, t):
if t < 2:
return model.Gaswarm[t] == 0
elif (model.Gaswarm[t-2] + model.Gaswarm[t-1] + model.Gascold[t-1]) == 1:
return model.Gaswarm[t] >= 0
else:
return model.Gaswarm[t] == 0
model.gas_warm_check = Constraint(model.T, rule=gas_warm)
"gas can only get hot from warm"
def gas_hot(model, t):
if t <= 2:
return model.Gashot[t] == 0
elif value(model.Gaswarm[t-2] + model.Gaswarm[t-1] + model.Gashot[t-2] + model.Gashot[t-1]) < 2:
return model.Gashot[t] == 0
else:
return model.Gashot[t] >= 0
model.gas_hot_check = Constraint(model.T, rule=gas_hot)
"Gas hot min load"
def gas_hot_min_load(model, t):
if model.Gashot[t] == 1:
return model.GasUnit_mw[t] >= 30 / 12
else:
return model.GasUnit_mw[t] == 0
model.gas_hot_min_load_check = Constraint(model.T, rule=gas_hot_min_load)
"Gas volume used"
def gas_volume_used(model, t):
'how much gas is being used?'
return model.GasUse[t] == (model.GasUnit_mw[t] * hr * (1+sys_use)) + (model.Gaswarm[t] * gunits * cold_gas_p_unit * (1+sys_use))
model.gas_vol = Constraint(model.T, rule=gas_volume_used)
"only 9000 GJ of gas can be used in a single day"
def daily_max_gas(model, t):
min_t = model.T.first()
max_t = model.T.last()
# Check all t until the last 24 hours
if t <= max_t:
return sum(model.GasUse[i] for i in range(min_t, max_t+1)) <= model.MaxVol
else:
return Constraint.Skip
model.max_gas_vol = Constraint(model.T, rule=daily_max_gas)
# Define the income, expenses, and profit
gas_income = sum(df.loc[t, 'rrp'] * model.GasUnit_mw[t] for t in model.T)
gas_expenses = sum(model.GasUse[t] * gprice + model.GasUnit_mw[t] * vom for t in model.T)
profit = gas_income - gas_expenses
model.objective = Objective(expr=profit, sense=maximize)
# Solve the model
solver = SolverFactory('glpk')
solver.solve(model)
results_df = model_to_df(model, first_period=first_model_period, last_period=last_model_period)
results_df['local_time'] = df.loc[:, 'local_time']
return results_df
# *********************
# RUN SCRIPT
# *********************
# Import file
df = pd.read_csv('raw_prices.csv')
# New Index
df['period'] = df.reset_index().index
# get current Year
first_model_period = 0
last_model_period = df['period'].iloc[-1]
# Gas Inputs
gunits = 3 # No units
gnet = 58.5 # MW/unit
srmc = 115.49 # $/MWh
vom = 6.18 # $/MWh
gprice = 10.206 # $/GJ
cold_gas_p_unit = 15.75 # GJ/unit/5mins
hr = 10.5 # GJ/MWh
sys_use = 0.02 # %
max_gas = 9000 # GJ
# Run Model
results_df = optimize_year(df=df, gunits=gunits, gnet=gnet, gprice=gprice, vom=vom, cold_gas_p_unit=cold_gas_p_unit,
hr=hr, max_gas=max_gas, sys_use=sys_use,
first_model_period=first_model_period, last_model_period=last_model_period)
# ---------------------
# Write the file out and remove the index
# ---------------------
print('ran successfully')
results_df.to_csv('optimised_output.csv', index=False)
示例数据
local_time rrp
3/07/2022 16:40 64.99
3/07/2022 16:45 73.10744
3/07/2022 16:50 74.21746
3/07/2022 16:55 95.90964
3/07/2022 17:00 89.60728
3/07/2022 17:05 90.37018
3/07/2022 17:10 92.52947
3/07/2022 17:15 100.1449
3/07/2022 17:20 103.70199
3/07/2022 17:25 290.2848
3/07/2022 17:30 115.40232
3/07/2022 17:35 121.00009
3/07/2022 17:40 294.62543
3/07/2022 17:45 121
3/07/2022 17:50 294.69446
3/07/2022 17:55 124.90664
3/07/2022 18:00 117.21929
3/07/2022 18:05 115.60193
3/07/2022 18:10 121
3/07/2022 18:15 148.20186
3/07/2022 18:20 123.11763
3/07/2022 18:25 120.99996
3/07/2022 18:30 121
3/07/2022 18:35 121
3/07/2022 18:40 121
3/07/2022 18:45 120.77726
3/07/2022 18:50 105.44435
3/07/2022 18:55 104.22928
3/07/2022 19:00 102.26646
3/07/2022 19:05 99.5896
3/07/2022 19:10 95.5881
3/07/2022 19:15 88
3/07/2022 19:20 88
3/07/2022 19:25 88
3/07/2022 19:30 96.97537
3/07/2022 19:35 88
3/07/2022 19:40 86.39322
3/07/2022 19:45 88
3/07/2022 19:50 85.97703
3/07/2022 19:55 86
3/07/2022 20:00 87.61047
3/07/2022 20:05 88
3/07/2022 20:10 88
3/07/2022 20:15 85.34908
3/07/2022 20:20 85.1
3/07/2022 20:25 71.01181
3/07/2022 20:30 84.88422
3/07/2022 20:35 85.39669
3/07/2022 20:40 86.02542
3/07/2022 20:45 84.87425
3/07/2022 20:50 70.0652
3/07/2022 20:55 70.65106
3/07/2022 21:00 70.53063
3/07/2022 21:05 68.27
3/07/2022 21:10 68.86159
3/07/2022 21:15 68.32018
3/07/2022 21:20 67.30473
3/07/2022 21:25 66.78292
3/07/2022 21:30 66.02195
3/07/2022 21:35 66.40755
好吧，我对此有点怪异。上钩了，挺有意思的问题。
因此，我进行了一系列更改，并在此示例中保留了您的一些代码。我还砍掉了一些成本变量并使它们变得相当简单，因为我有点迷失了酱汁，所以我（大部分）确信一切正常，所以 units/conversions/costs 现在有点荒谬, 但应该很容易恢复。
希望这里有几个概念供您在完成此过程时使用。一些注意事项...
- 需要一个二进制变量来指示工厂已启动，以及另一个来跟踪它是否在特定时期“运行ning”，这些与约束相关联
- 在时间 windows 上添加了一些小技巧，以便对总的 gas 使用量进行滚动评估
- 将工厂的最低使用量添加到 运行，否则一旦它“启动”，它可以任意 运行 在不盈利时使用 0 gas，现在 minimum-run 或关闭被迫的决定
情节显示非常有说服力的证据表明它如所希望的那样 运行ning，它启动，运行s 在价格高时达到最大爆炸，并遵守滚动气体限制，然后关闭然后再做一次。
代码
import pandas as pd
import numpy as np
from pyomo.environ import *
#from pyomo.core import * # not needed
import matplotlib.pyplot as plt
# Import file
df = pd.read_csv('raw_prices_full.csv')
# New Index
df['period'] = df.reset_index().index
print(df.head())
print(df.info())
# get current Year
first_model_period = 0
last_model_period = df['period'].iloc[-1]
# Gas Inputs
gunits = 2 # No units
gnet = 20 # MW/unit
#srmc = 115.49 # $/MWh
#vom = 6.18 # $/MWh
gprice = 10 # $/GJ
startup_gas_rate = 8 # GJ/unit/5mins
min_running = 5 # GJ/unit/5mins <--- assumed minimal output... or plant could "run" at zero
hr = 10 # GJ/MWh
sys_use = 0.10 # % overhead in GJ of gas
max_gas = 100 # GJ/30 minutes (6 periods), cut down for testing
# Filter the data
df = df.loc[first_model_period:last_model_period, :]
model = ConcreteModel()
# Define model parameters
model.T = Set(initialize=df.period.tolist(), doc='periods of year', ordered=True)
# P needs to be indexed by T, and should be init by dictionary
model.P = Param(model.T, initialize=df.rrp.to_dict(), doc='price for each period', within=NonNegativeReals)
model.MaxGas = Param(initialize=max_gas, doc='Max gas usage (GJ) in 24 hour period', within=NonNegativeReals)
# variables
model.turn_on = Var(model.T, domain=Binary) # decision to turn the plant on in period t, making it available at t+2
model.plant_running = Var(model.T, domain=Binary) # plant is running (able to produce heat) in period t. init is not necessary
model.total_gas = Var(model.T, bounds=(0, gunits*gnet / 12 * hr * (1+sys_use))) # the total gas flow in any period [GJ/5min]
model.gas_for_power = Var(model.T, domain=NonNegativeReals ) # the gas flow for power conversion at any period [GJ/5min]
model.profit = Var(model.T) # non-essential variable, just used for bookeeping
# *********************
# GAS CONSTRAINTS
# *********************
# logic for the plant_running control variable
def running(model, t):
# plant cannot be running in first 2 periods
if t < 2:
return model.plant_running[t] == 0
else:
# plant can be running if it was running in last period or started 2 periods ago.
return model.plant_running[t] <= model.plant_running[t-1] + model.turn_on[t-2]
model.running_constraint = Constraint(model.T, rule=running)
# charge the warmup gas after a start decision. Note the conditions in this constraint are all
# mutually exclusive, so there shouldn't be any double counting
# this will constrain the minimum gas flow to either the startup flow for 2 periods,
# the minimum running flow, or zero if neither condition is met
def gas_mins(model, t):
# figure out the time periods to inspect to see if a start event occurred...
if t==0:
possible_starts = model.turn_on[t]
else:
possible_starts = model.turn_on[t] + model.turn_on[t-1]
return model.total_gas[t] >= \
startup_gas_rate * gunits * (1 + sys_use) * possible_starts + \
min_running * gunits * (1 + sys_use) * model.plant_running[t]
model.gas_mins = Constraint(model.T, rule=gas_mins)
# constrain the gas used for power to zero if not running, or the period max
def max_flow(model, t):
return model.gas_for_power[t] <= model.plant_running[t] * gunits * gnet / 12 * hr
model.running_max_gas = Constraint(model.T, rule=max_flow)
# add the overhead to running gas to capture total gas by when running
def tot_gas_for_conversion(model, t):
return model.total_gas[t] >= model.gas_for_power[t] * (1 + sys_use)
model.tot_gas_for_conversion = Constraint(model.T, rule=tot_gas_for_conversion)
# I don't believe the constraint below is needed. The model will try to maximize
# gas use when profit is available, and we have already established a minimum
# def gas_volume_used(model, t):
# 'how much gas is being used?'
# 'during cold start, a small amount of gas is used'
# 'during operation, gas is consumed based on the Heat Rate'
# if model.plant_running[t] == 1 and model.GasUnit_mw[t] == 0:
# return model.gas_use[t] == startup_gas_rate * gunits * (1+sys_use)
# else:
# return model.gas_use[t] == model.GasUnit_mw[t] * hr * (1+sys_use)
# model.gas_vol = Constraint(model.T, rule=gas_volume_used)
# this isn't doing what you think, you are not controlling the range of summation to 24 time periods
# def daily_max_gas(model, t):
# "only 9000 GJ of gas can be used in a single day"
# min_t = model.T.first()
# max_t = model.T.last()
# # Check all t until the last 24 hours
# if t <= max_t:
# return sum(model.gas_use[i] for i in range(min_t, max_t+1)) <= model.MaxGas
# else:
# return Constraint.Skip
# model.max_gas_vol = Constraint(model.T, rule=daily_max_gas)
# constraint to limit consumption to max per 24hr rolling period. Note, this is *rolling* constraint,
# if a "daily" constraint tied to particular time of day is needed, more work will need to be done
# on the index to identify the indices of interest
window_size = 6 # for testing on small data, normally would be: 24hrs * 12intervals/hr
def rolling_max(model, t):
preceding_periods = {t_prime for t_prime in model.T if t - window_size <= t_prime < t}
return sum(model.total_gas[t_prime] for t_prime in preceding_periods) <= model.MaxGas
eval_periods = {t for t in model.T if t >= window_size}
model.rolling_max = Constraint(eval_periods, rule=rolling_max)
# Define the income, expenses, and profit
gas_income = sum(df.loc[t, 'rrp'] * model.gas_for_power[t] / hr for t in model.T)
gas_expenses = sum(model.total_gas[t] * gprice for t in model.T) # removed the "vom" computation
profit = gas_income - gas_expenses
model.objective = Objective(expr=profit, sense=maximize)
# Solve the model
solver = SolverFactory('glpk')
results = solver.solve(model)
print(results)
# model.display()
# model.pprint()
cols = []
cols.append(pd.Series(model.P.extract_values(), name='energy price'))
cols.append(pd.Series(model.total_gas.extract_values(), name='total gas'))
cols.append(pd.Series(model.gas_for_power.extract_values(), name='converted gas'))
df_results = pd.concat(cols, axis=1)
fig, ax1 = plt.subplots()
width = 0.25
ax1.bar(np.array(df.index)-width, df_results['energy price'], width, color='g', label='price of power')
ax1.set_ylabel('price')
ax2 = ax1.twinx()
ax2.set_ylabel('gas')
ax2.bar(df.index, df_results['total gas'], width, color='b', label='tot gas')
ax2.bar(np.array(df.index)+width, df_results['converted gas'], width, color='xkcd:orange', label='converted gas')
fig.legend()
fig.tight_layout()
plt.show()
背景
我正在尝试编写一个 pyomo 脚本，以根据对电价的完美预见来优化调度天然气厂。我相信我已经完成了 90%，只有几个问题。
问题
我的脚本有效，但求解器从不调度工厂，即使它应该在的地方，在下面提供的示例中，我可以手动计算至少 8131 美元的潜在利润。
我怀疑结果为零的原因是我编写约束的方式造成的，其中有 2 个；
- Gas Plant 从冷启动开始需要 10 分钟
- 预热后，燃气设备有一个必须运行的最小负载 at/above。
- 天然气厂一天只能消耗 9000 GJ 的天然气
特别是在进一步测试中，我认为是 'gas_volume_used' 约束导致了问题。
请求帮助
有人可以看看我的代码，看看我在约束方程中遗漏了什么吗？
代码
import pandas as pd
from pyomo.environ import *
from pyomo.core import *
# ---------------------
# Main Functions
# ---------------------
def model_to_df(model, first_period, last_period):
# Need to increase the first & last period by 1 because of pyomo indexing
# and add 1 to the value of last model period because of the range
periods = range(model.T[first_period + 1], model.T[last_period + 1] + 1)
rrp = [model.P.extract_values()[None][i] for i in periods]
Gascold = [value(model.Gascold[i]) for i in periods]
Gaswarm = [value(model.Gaswarm[i]) for i in periods]
Gashot = [value(model.Gashot[i]) for i in periods]
GasUnit_mw = [value(model.GasUnit_mw[i]) for i in periods]
GasUse = [value(model.GasUse[i]) for i in periods]
df_dict = dict(
period=periods,
rrp=rrp,
Gascold=Gascold,
Gaswarm=Gaswarm,
Gashot=Gashot,
GasUnit_mw=GasUnit_mw,
GasUse=GasUse
)
df = pd.DataFrame(df_dict)
return df
def optimize_year(df, gunits, gnet, gprice, vom, cold_gas_p_unit, hr, max_gas, sys_use,
first_model_period, last_model_period):
"""
Optimize behavior of a gas plant over a time period
Assume perfect foresight of electricity prices.
Parameters
----------
df : dataframe
dataframe with columns of hourly LBMP and the hour of the year
gunits : int
number of gas turbines
gnet : double
size in MW's of each gas turbine
gprice : double
the price of gas $/GJ
vom : double
the variable cost of gas $/MWh dispatched
cold_gas_p_unit : double
the amount of gas in GJ used in 5 minutes while booting up from cold start
hr : double
heat rate, the amount of gas consumed per MWh produced (GJ/MWh)
sys_use : double
auxillary gas usage, percentage
max_gas : double
the maximum amount of gas (GJ) available to use in a period
first_model_period : int, optional
Set the first period of the year to be considered in the optimization
last_model_period : int, optional
Set the last period of the year to be considered in the optimization
Returns
-------
dataframe
Gas plant operational stats and time stamp
"""
# Filter the data
df = df.loc[first_model_period:last_model_period, :]
model = ConcreteModel()
# Define model parameters
model.T = Set(initialize=df.period.tolist(), doc='periods of year', ordered=True)
model.P = Param(initialize=df.rrp.tolist(), doc='price for each period', within=Any)
model.MaxGas = Param(initialize=(gnet/12) * gunits * hr * (1+sys_use), doc='Max gas usage (GJ) in interval', within=Any)
model.MaxVol = Param(initialize=max_gas, doc='Max gas usage (GJ) in 24 hour period', within=Any)
model.MaxMW = Param(initialize=(gnet/12) * gunits, doc='Max MWs', within=Any)
model.Gascold = Var(model.T, doc='binary', within=Binary, initialize=1)
model.Gaswarm = Var(model.T, doc='binary', within=Binary, initialize=0)
model.Gashot = Var(model.T, doc='binary', within=Binary, initialize=0)
model.GasUse = Var(model.T, bounds=(0, model.MaxGas))
model.GasUnit_mw = Var(model.T, bounds=(0, model.MaxMW))
# *********************
# GAS CONSTRAINTS
# *********************
"Ensure that all three Gas State flags sum to 1"
def set_on(model, t):
return model.Gascold[t] + model.Gaswarm[t] + model.Gashot[t] == 1
model.limit_flags = Constraint(model.T, rule=set_on)
"Set t0 to GAS COLD = 1"
def gas_cold(model, t):
if t == model.T.first():
return model.Gascold[t] == 1
else:
return model.Gascold[t] >= 0
model.gas_cold_check = Constraint(model.T, rule=gas_cold)
"Gas can only warm up from cold"
def gas_warm(model, t):
if t < 2:
return model.Gaswarm[t] == 0
elif (model.Gaswarm[t-2] + model.Gaswarm[t-1] + model.Gascold[t-1]) == 1:
return model.Gaswarm[t] >= 0
else:
return model.Gaswarm[t] == 0
model.gas_warm_check = Constraint(model.T, rule=gas_warm)
"gas can only get hot from warm"
def gas_hot(model, t):
if t <= 2:
return model.Gashot[t] == 0
elif value(model.Gaswarm[t-2] + model.Gaswarm[t-1] + model.Gashot[t-2] + model.Gashot[t-1]) < 2:
return model.Gashot[t] == 0
else:
return model.Gashot[t] >= 0
model.gas_hot_check = Constraint(model.T, rule=gas_hot)
"Gas hot min load"
def gas_hot_min_load(model, t):
if model.Gashot[t] == 1:
return model.GasUnit_mw[t] >= 30 / 12
else:
return model.GasUnit_mw[t] == 0
model.gas_hot_min_load_check = Constraint(model.T, rule=gas_hot_min_load)
"Gas volume used"
def gas_volume_used(model, t):
'how much gas is being used?'
return model.GasUse[t] == (model.GasUnit_mw[t] * hr * (1+sys_use)) + (model.Gaswarm[t] * gunits * cold_gas_p_unit * (1+sys_use))
model.gas_vol = Constraint(model.T, rule=gas_volume_used)
"only 9000 GJ of gas can be used in a single day"
def daily_max_gas(model, t):
min_t = model.T.first()
max_t = model.T.last()
# Check all t until the last 24 hours
if t <= max_t:
return sum(model.GasUse[i] for i in range(min_t, max_t+1)) <= model.MaxVol
else:
return Constraint.Skip
model.max_gas_vol = Constraint(model.T, rule=daily_max_gas)
# Define the income, expenses, and profit
gas_income = sum(df.loc[t, 'rrp'] * model.GasUnit_mw[t] for t in model.T)
gas_expenses = sum(model.GasUse[t] * gprice + model.GasUnit_mw[t] * vom for t in model.T)
profit = gas_income - gas_expenses
model.objective = Objective(expr=profit, sense=maximize)
# Solve the model
solver = SolverFactory('glpk')
solver.solve(model)
results_df = model_to_df(model, first_period=first_model_period, last_period=last_model_period)
results_df['local_time'] = df.loc[:, 'local_time']
return results_df
# *********************
# RUN SCRIPT
# *********************
# Import file
df = pd.read_csv('raw_prices.csv')
# New Index
df['period'] = df.reset_index().index
# get current Year
first_model_period = 0
last_model_period = df['period'].iloc[-1]
# Gas Inputs
gunits = 3 # No units
gnet = 58.5 # MW/unit
srmc = 115.49 # $/MWh
vom = 6.18 # $/MWh
gprice = 10.206 # $/GJ
cold_gas_p_unit = 15.75 # GJ/unit/5mins
hr = 10.5 # GJ/MWh
sys_use = 0.02 # %
max_gas = 9000 # GJ
# Run Model
results_df = optimize_year(df=df, gunits=gunits, gnet=gnet, gprice=gprice, vom=vom, cold_gas_p_unit=cold_gas_p_unit,
hr=hr, max_gas=max_gas, sys_use=sys_use,
first_model_period=first_model_period, last_model_period=last_model_period)
# ---------------------
# Write the file out and remove the index
# ---------------------
print('ran successfully')
results_df.to_csv('optimised_output.csv', index=False)
示例数据
local_time rrp
3/07/2022 16:40 64.99
3/07/2022 16:45 73.10744
3/07/2022 16:50 74.21746
3/07/2022 16:55 95.90964
3/07/2022 17:00 89.60728
3/07/2022 17:05 90.37018
3/07/2022 17:10 92.52947
3/07/2022 17:15 100.1449
3/07/2022 17:20 103.70199
3/07/2022 17:25 290.2848
3/07/2022 17:30 115.40232
3/07/2022 17:35 121.00009
3/07/2022 17:40 294.62543
3/07/2022 17:45 121
3/07/2022 17:50 294.69446
3/07/2022 17:55 124.90664
3/07/2022 18:00 117.21929
3/07/2022 18:05 115.60193
3/07/2022 18:10 121
3/07/2022 18:15 148.20186
3/07/2022 18:20 123.11763
3/07/2022 18:25 120.99996
3/07/2022 18:30 121
3/07/2022 18:35 121
3/07/2022 18:40 121
3/07/2022 18:45 120.77726
3/07/2022 18:50 105.44435
3/07/2022 18:55 104.22928
3/07/2022 19:00 102.26646
3/07/2022 19:05 99.5896
3/07/2022 19:10 95.5881
3/07/2022 19:15 88
3/07/2022 19:20 88
3/07/2022 19:25 88
3/07/2022 19:30 96.97537
3/07/2022 19:35 88
3/07/2022 19:40 86.39322
3/07/2022 19:45 88
3/07/2022 19:50 85.97703
3/07/2022 19:55 86
3/07/2022 20:00 87.61047
3/07/2022 20:05 88
3/07/2022 20:10 88
3/07/2022 20:15 85.34908
3/07/2022 20:20 85.1
3/07/2022 20:25 71.01181
3/07/2022 20:30 84.88422
3/07/2022 20:35 85.39669
3/07/2022 20:40 86.02542
3/07/2022 20:45 84.87425
3/07/2022 20:50 70.0652
3/07/2022 20:55 70.65106
3/07/2022 21:00 70.53063
3/07/2022 21:05 68.27
3/07/2022 21:10 68.86159
3/07/2022 21:15 68.32018
3/07/2022 21:20 67.30473
3/07/2022 21:25 66.78292
3/07/2022 21:30 66.02195
3/07/2022 21:35 66.40755
好吧，我对此有点怪异。上钩了，挺有意思的问题。
因此，我进行了一系列更改，并在此示例中保留了您的一些代码。我还砍掉了一些成本变量并使它们变得相当简单，因为我有点迷失了酱汁，所以我（大部分）确信一切正常，所以 units/conversions/costs 现在有点荒谬, 但应该很容易恢复。
希望这里有几个概念供您在完成此过程时使用。一些注意事项...
- 需要一个二进制变量来指示工厂已启动，以及另一个来跟踪它是否在特定时期“运行ning”，这些与约束相关联
- 在时间 windows 上添加了一些小技巧，以便对总的 gas 使用量进行滚动评估
- 将工厂的最低使用量添加到 运行，否则一旦它“启动”，它可以任意 运行 在不盈利时使用 0 gas，现在 minimum-run 或关闭被迫的决定
情节显示非常有说服力的证据表明它如所希望的那样 运行ning，它启动，运行s 在价格高时达到最大爆炸，并遵守滚动气体限制，然后关闭然后再做一次。
代码
import pandas as pd
import numpy as np
from pyomo.environ import *
#from pyomo.core import * # not needed
import matplotlib.pyplot as plt
# Import file
df = pd.read_csv('raw_prices_full.csv')
# New Index
df['period'] = df.reset_index().index
print(df.head())
print(df.info())
# get current Year
first_model_period = 0
last_model_period = df['period'].iloc[-1]
# Gas Inputs
gunits = 2 # No units
gnet = 20 # MW/unit
#srmc = 115.49 # $/MWh
#vom = 6.18 # $/MWh
gprice = 10 # $/GJ
startup_gas_rate = 8 # GJ/unit/5mins
min_running = 5 # GJ/unit/5mins <--- assumed minimal output... or plant could "run" at zero
hr = 10 # GJ/MWh
sys_use = 0.10 # % overhead in GJ of gas
max_gas = 100 # GJ/30 minutes (6 periods), cut down for testing
# Filter the data
df = df.loc[first_model_period:last_model_period, :]
model = ConcreteModel()
# Define model parameters
model.T = Set(initialize=df.period.tolist(), doc='periods of year', ordered=True)
# P needs to be indexed by T, and should be init by dictionary
model.P = Param(model.T, initialize=df.rrp.to_dict(), doc='price for each period', within=NonNegativeReals)
model.MaxGas = Param(initialize=max_gas, doc='Max gas usage (GJ) in 24 hour period', within=NonNegativeReals)
# variables
model.turn_on = Var(model.T, domain=Binary) # decision to turn the plant on in period t, making it available at t+2
model.plant_running = Var(model.T, domain=Binary) # plant is running (able to produce heat) in period t. init is not necessary
model.total_gas = Var(model.T, bounds=(0, gunits*gnet / 12 * hr * (1+sys_use))) # the total gas flow in any period [GJ/5min]
model.gas_for_power = Var(model.T, domain=NonNegativeReals ) # the gas flow for power conversion at any period [GJ/5min]
model.profit = Var(model.T) # non-essential variable, just used for bookeeping
# *********************
# GAS CONSTRAINTS
# *********************
# logic for the plant_running control variable
def running(model, t):
# plant cannot be running in first 2 periods
if t < 2:
return model.plant_running[t] == 0
else:
# plant can be running if it was running in last period or started 2 periods ago.
return model.plant_running[t] <= model.plant_running[t-1] + model.turn_on[t-2]
model.running_constraint = Constraint(model.T, rule=running)
# charge the warmup gas after a start decision. Note the conditions in this constraint are all
# mutually exclusive, so there shouldn't be any double counting
# this will constrain the minimum gas flow to either the startup flow for 2 periods,
# the minimum running flow, or zero if neither condition is met
def gas_mins(model, t):
# figure out the time periods to inspect to see if a start event occurred...
if t==0:
possible_starts = model.turn_on[t]
else:
possible_starts = model.turn_on[t] + model.turn_on[t-1]
return model.total_gas[t] >= \
startup_gas_rate * gunits * (1 + sys_use) * possible_starts + \
min_running * gunits * (1 + sys_use) * model.plant_running[t]
model.gas_mins = Constraint(model.T, rule=gas_mins)
# constrain the gas used for power to zero if not running, or the period max
def max_flow(model, t):
return model.gas_for_power[t] <= model.plant_running[t] * gunits * gnet / 12 * hr
model.running_max_gas = Constraint(model.T, rule=max_flow)
# add the overhead to running gas to capture total gas by when running
def tot_gas_for_conversion(model, t):
return model.total_gas[t] >= model.gas_for_power[t] * (1 + sys_use)
model.tot_gas_for_conversion = Constraint(model.T, rule=tot_gas_for_conversion)
# I don't believe the constraint below is needed. The model will try to maximize
# gas use when profit is available, and we have already established a minimum
# def gas_volume_used(model, t):
# 'how much gas is being used?'
# 'during cold start, a small amount of gas is used'
# 'during operation, gas is consumed based on the Heat Rate'
# if model.plant_running[t] == 1 and model.GasUnit_mw[t] == 0:
# return model.gas_use[t] == startup_gas_rate * gunits * (1+sys_use)
# else:
# return model.gas_use[t] == model.GasUnit_mw[t] * hr * (1+sys_use)
# model.gas_vol = Constraint(model.T, rule=gas_volume_used)
# this isn't doing what you think, you are not controlling the range of summation to 24 time periods
# def daily_max_gas(model, t):
# "only 9000 GJ of gas can be used in a single day"
# min_t = model.T.first()
# max_t = model.T.last()
# # Check all t until the last 24 hours
# if t <= max_t:
# return sum(model.gas_use[i] for i in range(min_t, max_t+1)) <= model.MaxGas
# else:
# return Constraint.Skip
# model.max_gas_vol = Constraint(model.T, rule=daily_max_gas)
# constraint to limit consumption to max per 24hr rolling period. Note, this is *rolling* constraint,
# if a "daily" constraint tied to particular time of day is needed, more work will need to be done
# on the index to identify the indices of interest
window_size = 6 # for testing on small data, normally would be: 24hrs * 12intervals/hr
def rolling_max(model, t):
preceding_periods = {t_prime for t_prime in model.T if t - window_size <= t_prime < t}
return sum(model.total_gas[t_prime] for t_prime in preceding_periods) <= model.MaxGas
eval_periods = {t for t in model.T if t >= window_size}
model.rolling_max = Constraint(eval_periods, rule=rolling_max)
# Define the income, expenses, and profit
gas_income = sum(df.loc[t, 'rrp'] * model.gas_for_power[t] / hr for t in model.T)
gas_expenses = sum(model.total_gas[t] * gprice for t in model.T) # removed the "vom" computation
profit = gas_income - gas_expenses
model.objective = Objective(expr=profit, sense=maximize)
# Solve the model
solver = SolverFactory('glpk')
results = solver.solve(model)
print(results)
# model.display()
# model.pprint()
cols = []
cols.append(pd.Series(model.P.extract_values(), name='energy price'))
cols.append(pd.Series(model.total_gas.extract_values(), name='total gas'))
cols.append(pd.Series(model.gas_for_power.extract_values(), name='converted gas'))
df_results = pd.concat(cols, axis=1)
fig, ax1 = plt.subplots()
width = 0.25
ax1.bar(np.array(df.index)-width, df_results['energy price'], width, color='g', label='price of power')
ax1.set_ylabel('price')
ax2 = ax1.twinx()
ax2.set_ylabel('gas')
ax2.bar(df.index, df_results['total gas'], width, color='b', label='tot gas')
ax2.bar(np.array(df.index)+width, df_results['converted gas'], width, color='xkcd:orange', label='converted gas')
fig.legend()
fig.tight_layout()
plt.show()