Pyomo 优化不起作用(天然气厂调度)

Pyomo optimisation not working (gas plant dispatch)

背景

我正在尝试编写一个 pyomo 脚本,以根据对电价的完美预见来优化调度天然气厂。我相信我已经完成了 90%,只有几个问题。

问题

我的脚本有效,但求解器从不调度工厂,即使它应该在的地方,在下面提供的示例中,我可以手动计算至少 8131 美元的潜在利润。

我怀疑结果为零的原因是我编写约束的方式造成的,其中有 2 个;

  1. Gas Plant 从冷启动开始需要 10 分钟
  2. 预热后,燃气设备有一个必须运行的最小负载 at/above。
  3. 天然气厂一天只能消耗 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()