Pyomo 优化不起作用（天然气厂调度）

Question

背景

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

问题

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

我怀疑我的零结果的原因是由于我如何编写约束，其中有 2 个；

1. 燃气厂需要 10 分钟才能从冷启动启动。
2. 一旦预热，燃气厂就有一个必须在/以上运行的最小负载。
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

Answer 1

好吧，我对此有点极客。上瘾了，有点有趣的问题。

因此，我做了很多更改，并在此示例中留下了一些代码。我还削减了一些成本变量并使它们变得相当简单，因为我有点迷失在酱汁中，所以我（大部分）确信事情正在发挥作用，所以单位/转换/成本现在有点荒谬，但应该很容易恢复。

希望这里有一些概念可供您在完成此操作时使用。一些注释...

• 需要一个二进制变量来指示工厂已启动，另一个变量来跟踪它是否在特定时期“运行”，这些变量与约束相关联
• 添加了一些与时间有关的技巧为总天然气使用量进行滚动评估期的窗口
• 添加了工厂运行的最低使用量，否则一旦“启动”，它可以在无利可图时任意使用 0 天然气运行，现在强制执行最低运行或关闭

决定显示了非常令人信服的证据表明它正在运行正如所希望的那样，它启动，在价格高时以最大速度运行，并遵守滚动燃气限制，然后关闭并再次执行。

代码

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()