Pyomo中如何构造二维Set的子集以简化BEV充电优化约束

我需要求解沿电网线路的纯电动汽车(BEV)最优充电策略，因此使用pyomo搭建描述该优化问题的线性规划(LP)模型。我首先定义了步长为15分钟的全天时间集合：

model.times = pe.Set(initialize=list(range(96)))

以及电网节点集合：

model.buses = pe.Set(initialize=list(range(6)))

我使用这两个集合的笛卡尔积索引变量和约束，例如跟踪BEV荷电状态(SOC)的约束如下：

def track_socs_rule(model, t, b):
    if t >= model.bevs[b].t_start and t <= model.bevs[b].t_target:
        return (model.SOC[t, b] + model.I[t, b] * model.voltages[b] * model.resolution/60 / 1000
                / model.bevs[b].e_bat*100 - model.SOC[t+1, b]) == 0
    else:
        return pe.Constraint.Skip

当前我需要在约束定义规则中手动通过条件判断过滤有效时间步（即对应节点的BEV处于充电时段的时间步），虽然该实现可得到正确结果，但代码写法不够优雅。我尝试直接定义仅包含有效(时间,节点)对的二维子集，写法如下：

model.valid_subset = pe.Set(within=model.times*model.charger_buses,
                            initialize=[[t for t in model.times if t >= model.bevs[b].t_start and
                                         t <= model.bevs[b].t_target] for b in model.charger_buses]

但运行时抛出如下错误：

ValueError: Cannot add value([bla, bla, ...], [bla, bla...], ...) to Set valid_subset.
The value is not in the domain valid_subset_domain

以下是可复现该错误的最小示例代码：

import pyomo.environ as pe

class BEV:
    def __init__(self, home_bus, e_bat, soc_start, soc_target, t_start, t_target, resolution):
        self.resolution = resolution
        self.home_bus = home_bus
        self.e_bat = e_bat
        self.soc_start = soc_start
        self.soc_target = soc_target
        self.t_start = int(t_start * 60/self.resolution)
        self.t_target = int(t_target * 60/self.resolution)


RESOLUTION = 15 # minutes

home_buses = [0, 2, 3, 4]
e_bats = [50, 50, 50, 50] # kWh
soc_starts = [20, 25, 20, 30] # %
soc_targets = [80, 75, 90, 100] # %
t_starts = [10, 12, 9, 15] # hours
t_targets = [19, 18, 15, 23] # hours

bevs = []
for pos in range(4):
    bev = BEV(home_bus=home_buses[pos], e_bat=e_bats[pos], soc_start=soc_starts[pos],
              soc_target=soc_targets[pos], t_start=t_starts[pos],
              t_target=t_targets[pos], resolution=RESOLUTION)
    bevs.append(bev)

solver = pe.SolverFactory('glpk')
model = pe.ConcreteModel()

model.bevs = {bev.home_bus: bev for bev in bevs}
model.resolution = RESOLUTION
model.times = pe.Set(initialize=list(range(int(24 * 60/RESOLUTION))))
model.buses = pe.Set(initialize=list(range(6)))
model.charger_buses = pe.Set(within=model.buses, initialize=[bus for bus in model.bevs])
#model.valid_subset = pe.Set(within=model.times*model.charger_buses,
                            #initialize=[[t for t in model.times if t >= model.bevs[b].t_start
                                         #and t <= model.bevs[b].t_target] for b in model.charger_buses])

model.voltages = pe.Param(model.buses, initialize={i: 400-i/2 for i in model.buses})

model.SOC = pe.Var(model.times*model.charger_buses, domain=pe.PositiveReals)
model.I = pe.Var(model.times*model.charger_buses, domain=pe.PositiveReals)

def track_socs_rule(model, t, b):
    # instead leave if (once indexed in model.valid_subset)
    if t >= model.bevs[b].t_start and t <= model.bevs[b].t_target:
        return (model.SOC[t, b] + model.I[t, b] * model.voltages[b] * model.resolution/60 /1000
                /model.bevs[b].e_bat*100 - model.SOC[t+1, b]) == 0

    else:
        return pe.Constraint.Skip

# instead index in model.valid_subset (once it works...)
model.track_socs = pe.Constraint(model.times*model.charger_buses, rule=track_socs_rule)
model.track_socs.pprint()

取消model.valid_subset相关代码的注释即可复现错误。我希望构造该二维有效子集后，可以直接用该子集索引约束，无需在约束规则中编写条件判断，同时也可用于变量I和SOC的上下界定义。请问在pyomo中是否存在正确构造二维Set子集的方法？

报错核心原因是传给initialize的是嵌套列表结构，每个元素是单个节点对应的有效时间列表，而二维Set要求传入的元素是(时间, 节点)格式的元组。

将valid_subset的定义替换为以下写法即可：

model.valid_subset = pe.Set(
    within=model.times * model.charger_buses,
    initialize= [
        (t, b) for b in model.charger_buses 
        for t in model.times 
        if t >= model.bevs[b].t_start and t <= model.bevs[b].t_target
    ]
)

构造完成后，可直接用该子集索引约束，简化后的约束定义如下：

def track_socs_rule(model, t, b):
    return (model.SOC[t, b] + model.I[t, b] * model.voltages[b] * model.resolution/60 /1000
            /model.bevs[b].e_bat*100 - model.SOC[t+1, b]) == 0

model.track_socs = pe.Constraint(model.valid_subset, rule=track_socs_rule)

无需再在规则内部写条件判断跳过非有效索引。

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