调试Pyomo中部分能量请求未被正常满足的问题

问题概述 #

基于Pyomo实现的能量分配模型，目标是在随时间变化的可用能量中最大化利用率，分配给带功率、能量需求和时间窗约束的请求。但出现反常现象：

• 单个请求需连续10个30分钟时段（共5小时）提供3.25kW，总能量16.25kWh，时间窗为910177-910211（35个时段）时，请求未被满足（m.satisfied == 0）；
• 将时间窗缩短为910177-910205（29个时段）时，请求反而能被满足（m.satisfied == 1），与“时间窗越小求解越难”的理论相悖。

核心问题定位 #

模型约束逻辑存在漏洞，未强制请求启动后连续运行指定时长：

• 原约束link_running仅限制m.running[t,r] <= m.running[t-1,r] + m.start[t,r]，允许启动后中途停止运行；
• 当时间窗较大时，模型可能生成断断续续的运行方案，总能量无法达到请求需求，导致satisfied=0；而时间窗较小时，仅存连续运行的可行解，因此能满足请求；
• 同时缺少对启动时间的边界约束：未确保启动后能在请求最晚结束时间前完成整个运行周期。

调试建议 #

1. 变量取值验证：在求解后输出m.start、m.running的具体取值，查看大时间窗下请求的运行是否为非连续状态；
2. 约束逻辑检查：验证每个约束是否准确匹配业务规则（如“必须连续运行指定时长”）；
3. 边界条件测试：手动计算请求的合法启动时间范围（启动时间t需满足t + duration -1 <= latest_end_time），检查模型是否允许超出该范围的启动。

优化后的模型实现 #

针对约束漏洞进行修改，补充连续运行和时间边界约束：

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 # 30分钟时间分辨率

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

# 可用功率供给（设为极大值排除供给限制影响）
supply_index = list(range(910176, 910223+1))
supply_vals = np.ones(len(supply_index))*2500000
supply_unix_ts = pd.Series(supply_vals, index=supply_index) 

# 请求参数
request_earliest_start_time = 910177
request_latest_end_time = 910211 # 修改为910205可复现原正常情况
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 # 需连续运行的时段数

# 构建每个时段的合格请求列表
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)

### 模型构建
m = pyo.ConcreteModel('dispatch')

### 集合
m.T = pyo.Set(initialize=tuple(range(910176,910223+1)))
m.R = pyo.Set(initialize=tuple(request_max_power_dict.keys()))
m.windows = pyo.Set(m.T, initialize=eligible_requests, within=m.R, doc='各时段合格请求')
m.windows_flat = pyo.Set(initialize={(t, r) for t in eligible_requests for r in eligible_requests[t]}, within=m.T * m.R)

### 参数
m.power_limit = pyo.Param(m.T, initialize=supply_unix_ts.to_dict()) # 各时段可用功率上限
m.request_power_size = pyo.Param(m.R, initialize=request_max_power_dict) # 请求功率需求
m.request_energy_size = pyo.Param(m.R, initialize=request_max_energy_dict) # 请求总能量需求
m.duration_periods = pyo.Param(m.R, initialize=request_duration_num_periods_dict) # 请求需连续运行的时段数
m.request_earliest_start = pyo.Param(m.R, initialize={r: s for r, (s, e) in request_time_constraints_dict.items()})
m.request_latest_end = pyo.Param(m.R, initialize={r: e for r, (s, e) in request_time_constraints_dict.items()})

### 变量
m.booked_supply = pyo.Var(m.T, domain=pyo.NonNegativeReals, doc='各时段已用功率')
m.dispatch = pyo.Var(m.windows_flat, domain=pyo.NonNegativeReals, doc='各时段分配给请求的功率')
m.start = pyo.Var(m.windows_flat, domain=pyo.Binary, doc='请求是否在该时段启动')
m.running = pyo.Var(m.windows_flat, domain=pyo.Binary, doc='请求是否在该时段运行')
m.satisfied = pyo.Var(m.R, domain=pyo.Binary, doc='请求是否被完全满足') 

### 目标函数：最大化能量利用率
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)

### 约束
# 1. 各时段已用功率不超过供给上限
@m.Constraint(m.T)
def supply_limit(m, t):
    return m.booked_supply[t] <= m.power_limit[t]

# 2. 各时段已用功率等于所有运行请求的功率之和
@m.Constraint(m.T)
def booked_supply_limit(m, t):
    return sum(m.running[t,r]*m.request_power_size[r] for r in m.windows[t]) == m.booked_supply[t]

# 3. 请求满足时，总能量等于需求；未满足时总能量为0
@m.Constraint(m.R)
def request_satisfied(m, r):
    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]

# 4. 每个请求最多启动一次
@m.Constraint(m.R)
def run_once(m, r):
    return sum(m.start[t,r] for t in m.T if r in m.windows[t]) <= 1

# 5. 强制请求启动后连续运行指定时长，且运行状态与启动状态关联
@m.Constraint(m.windows_flat)
def link_running(m, t, r):
    duration = int(m.duration_periods[r])
    # 启动时段：running必须等于start
    if (t, r) in m.windows_flat:
        # 若当前时段启动，则必须运行
        start_expr = m.start[t,r]
        # 检查当前时段是否在某个启动后的运行周期内
        in_run_window = sum(m.start[t0,r] for t0 in range(max(m.request_earliest_start[r], t - duration + 1), t) if (t0, r) in m.windows_flat and t0 + duration > t)
        return m.running[t,r] == start_expr + in_run_window - start_expr*in_run_window

# 6. 限制启动时间：启动后必须能在最晚结束时间前完成整个运行周期
@m.Constraint(m.windows_flat)
def valid_start_time(m, t, r):
    return m.start[t,r] <= (1 if t + m.duration_periods[r] -1 <= m.request_latest_end[r] else 0)

# 7. 跟踪实际分配功率（调试用）
@m.Constraint(m.windows_flat)
def keep_track_power(m, t, r):
    return m.running[t, r]*m.request_power_size[r] == m.dispatch[t, r]

# 求解模型
solver = pyo.SolverFactory('cbc')
res = solver.solve(m)

# 输出结果到文件
temp = sys.stdout 
f = open(str(int(time.time()))+'.txt', 'w')
sys.stdout = f
print("请求满足状态：", pyo.value(m.satisfied[0]))
print("启动时段：", [(t, pyo.value(m.start[t,0])) for t in m.T if (t,0) in m.windows_flat and pyo.value(m.start[t,0])>0])
print("运行时段：", [(t, pyo.value(m.running[t,0])) for t in m.T if (t,0) in m.windows_flat and pyo.value(m.running[t,0])>0])
m.pprint()
print(res)

f.close()
sys.stdout = temp

优化说明 #

1. 新增启动时间约束：valid_start_time确保请求启动后能在最晚结束时间前完成全部运行时段；
2. 修改运行状态约束：link_running强制请求启动后连续运行指定时长，禁止中途停止；
3. 新增参数：将请求的最早启动、最晚结束时间设为参数，简化约束逻辑；
4. 增强调试输出：直接打印请求的满足状态、启动和运行时段，方便验证结果。

内容的提问来源于stack exchange，提问作者user23390596