调试为什么有些请求在pyomo中意外得不到满足

问题描述 投票:0回答:1

我有一个能源分配问题,已在 Pyomo 中实施。目标是以最大化可用能量利用率的方式在一组请求之间分配可用能量(这是有限的且随时间变化的),该问题在此处进一步描述:https://stackoverflow.com/问题/77982332/生成约束与条件求和多个-pyomo-子集/77982847#77982847.

每个请求都有功率 (kW) 和能量 (kWh) 要求以及满足时间范围的限制。

我已经成功地实现了 99% 的请求(非常感谢@Airsquid),但我有几个案例,其中一些请求似乎从未得到满足,我无法弄清楚这是为什么。

我在下面复制了一个最小的例子。就上下文而言,时间集超过 24 小时,分辨率为 30 分钟(48 个周期,索引范围为 910176 到 910223(含))。

下例中的单个请求需要在 10 个周期(5 小时)内输送 3.25kW 的电力才能接收 16.25kWh,并且需要在 910177 到 910211(这是 35 个周期)的区间内以单个块接收该能量间隔或 17.5 小时)。

执行此代码时,请求未得到满足,即执行后 m.satisfied == 0。但是,如果我将最晚结束时间从 910211 更改为 910205(即,将传输能量的有效间隔从 35 个周期减少到 29 个周期),则请求确实会得到满足,即 m.satisfied == 1。

进行此更改对我来说没有意义,因为通过减少最晚结束时间,理论上我们会使问题更难解决。

如前所述,当我执行市场时,99% 的请求均按预期满足此代码,但即使当我执行仅包含有问题的请求的市场时,也有一些请求不满意。我无法确定与可能导致问题的请求有关的任何特定内容,例如,存在大/小功率请求、大/小能量请求、请求能量的大/小窗口的混合。

此时我非常不知道为什么会发生这种情况,我想知道这是否是我的实现中的问题,或者这是否可能是与 Pyomo 有关的问题。任何有关如何解决此问题的建议将不胜感激,也欢迎有关如何改进实施的最佳实践的一般建议!

import pyomo.environ as pyo
import sys
import time
from collections import defaultdict
import pandas as pd
import numpy as np

sampling_period_hrs=0.5 #our market uses a time resolution of 30 minutes / 0.5 hours, we use this later to relate request power delivered to request energy required

request_time_constraints_dict = {}
request_max_power_dict = {}
request_max_energy_dict={}
request_duration_num_periods_dict = {}

#this is the available supply (kW) to satisfy requests
supply_index = list(range(910176, 910223+1)) #this is the timeframe over which power is available to be allocated
supply_vals = np.ones(len(supply_index))*2500000 #this is the amount of power (kW) available, this has been made arbitrarily large for the purposes of this illustration
supply_unix_ts = pd.Series(supply_vals, index=supply_index) 


request_earliest_start_time = 910177 #earliest start time request can be satisfied
request_latest_end_time = 910211 #latest end time request can be satisfied, TEST: change this to 910205 to see request m.satisfied goes to 1 
request_time_constraints_dict[0] = (int(request_earliest_start_time), int(request_latest_end_time))
request_max_power_dict[0] = 3.25 #kw
request_max_energy_dict[0] = 16.25 #kwh
request_duration_num_periods_dict[0] = 10.0 #periods (10 x 0.5h periods = 5 hours)


eligible_requests = defaultdict(list)
for r, (start, end) in request_time_constraints_dict.items():
    
    for t in range(start, end + 1):
        eligible_requests[t].append(r)


### MODEL BUILD

m = pyo.ConcreteModel('dispatch')

### SETS
# m.T = pyo.Set(initialize=tuple(supply_unix_ts.index))
m.T = pyo.Set(initialize=tuple(range(910176,910223+1)))
m.R = pyo.Set(initialize=tuple(request_max_power_dict.keys()))
# Note:  This will be an "indexed" set, where we have sets indexed by some index, in this case, m.T
m.windows = pyo.Set(m.T, initialize=eligible_requests, within=m.R,
                    doc='requests eligible in each timeslot')
# We can make a flat-set here which will be a good basis for the dispatch variable
m.windows_flat = pyo.Set(initialize={(t, r) for t in eligible_requests for r in
                                    eligible_requests[t]},
                        within=m.T * m.R)


## PARAMS
m.power_limit = pyo.Param(m.T, initialize=supply_unix_ts.to_dict()) #total amount of power (kW) available for a specific period
m.request_power_size = pyo.Param(m.R, initialize=request_max_power_dict) #amount of power (kW) required for a request
m.request_energy_size = pyo.Param(m.R, initialize=request_max_energy_dict) #amount of energy (kWh) required for a request
m.duration_periods = pyo.Param(m.R, initialize=request_duration_num_periods_dict) #number of periods (note this is calculated by relating the energy to the power to ensure there are no conflicts)


### VARS
m.booked_supply = pyo.Var(m.T, domain = pyo.NonNegativeReals, doc='total amount of power used in a time period')
m.dispatch = pyo.Var(m.windows_flat, domain=pyo.NonNegativeReals,
                    doc='dispatch power in timeslot t to request r')
m.start = pyo.Var(m.windows_flat, domain=pyo.Binary, doc='binary for period in which request starts to be satisfied')
m.running = pyo.Var(m.windows_flat, domain=pyo.Binary, doc='binary for if request being satisfie in interval')
m.satisfied = pyo.Var(m.R, domain=pyo.Binary, doc='binary for if request is fully satisfied') 


### OBJ
def obj_rule(m):
        return sum(m.booked_supply[t] for t in m.T)
m.obj = pyo.Objective(rule = obj_rule, sense = pyo.maximize)


### CONSTRAINTS

#constraint 1
@m.Constraint(m.T)
def supply_limit(m, t):
    return m.booked_supply[t] <= m.power_limit[t]

#constraint 2
@m.Constraint(m.T)
def booked_supply_limit(m, t):
    # total power delivered in each interval
    return sum(m.running[t,r]*m.request_power_size[r] for r in m.windows[t]) == m.booked_supply[t]

#constraint 3
@m.Constraint(m.R)
def request_satisfied(m, r):
    # dictates whether a request is satisfied and relates power to energy
    timeslots = {t for t in m.T if r in m.windows[t]}
    return sum(m.running[t, r]*m.request_power_size[r] for t in timeslots)*sampling_period_hrs == m.satisfied[r] * m.request_energy_size[r]


# constraint 4
@m.Constraint(m.R)
def run_once(m, r):
    #request can only run once
    return(sum(m.start[t,r] for t in m.T if r in m.windows[t]) <= 1)

#constraint 5
@m.Constraint(m.windows_flat)
def link_running(m, t, r):
    #request can only be satisfied if it was running in prior period or started in this one (we need request to be satisfied in a single block rather than turning off/on)
    timeslots = {t for t in m.T if r in m.windows[t]}
    if t <= list(timeslots)[0]: 
        return m.running[t,r] == m.start[t,r]
    return m.running[t,r] <= m.running[t-1,r] + m.start[t,r]


# constraint 6
@m.Constraint(m.windows_flat)
def keep_track_power(m, t, r):
    #useful for debugging - keep track of actual power values for each request in each time interval
    return m.running[t, r]*m.request_power_size[r] == m.dispatch[t, r]


# solve it
solver = pyo.SolverFactory('cbc')
res = solver.solve(m)

temp = sys.stdout 
f = open(str(int(time.time()))+'.txt', 'w')
sys.stdout = f
m.pprint()
print(res)

f.close()
sys.stdout = temp
python debugging optimization mathematical-optimization pyomo
1个回答
0
投票

您的

link_running
约束中有一个错误,导致了所有的麻烦。如果您查看失败运行的一些输出,您会在该约束的打印输出中注意到,
running = start
的“第一周期”方程似乎到处都出现了......

那是因为您试图使用以下语句获取窗口中第一个符合条件的开始时间:

if t <= list(timeslots)[0]: 

timeslots
是一个 set,因此绝对不能保证 [0] 元素是最小值。你想这样做:

if t == min(timeslots):

正确识别第一个合格期间。


在第二个不相关的注释中,您可能可以重新考虑这一点并使其更简单在您使当前表单正常工作之后。如果需要的话。如果它运行良好且及时,那么就不用打扰。

但是,在您的模型中,似乎每个请求的电量在用于满足它的所有周期中都是“固定”的。真的吗?或者您是否打算创建一个更灵活的模型,其中每个请求分配的功率可以变化,因此用于满足请求的周期数可以变化?这是两种不同的模型。

在您当前的版本中,它似乎是前者。在这种情况下,您只需通过 start

变量就可以“了解有关请求的所有信息”。因此,您可以进行一次大的重构,并可能删除一些辅助变量。如果它从第 5 周期开始,并且请求持续 10 个周期且具有一定的固定电量,则您不需要变量来知道它是否在第 8 周期运行 -

它是!

 类似地,对于时期。您仍然需要一个变量来跟踪分配给每个请求的电量。
只是一些值得思考的事情--

© www.soinside.com 2019 - 2024. All rights reserved.