Pyomo 中的双索引变量和多个单索引变量有区别吗？

Question

我目前正在研究一个问题，我想在一段时间内优化多个对象（具体来说：几辆电动汽车在一段时间内的充电曲线）。

是否为每个电动汽车充电配置文件创建一个变量 (

pyo.Var

) 并使用时间步长对其进行索引，或者创建一个变量并使用时间步长和电动汽车对其进行索引，这有什么区别吗？或者换句话说：在扩大电动汽车数量或时间步数时，这种差异是否超出了个人品味，在性能方面是否有更好/更实用的变体？

我将尝试在 MWE 的帮助下说明差异：

底层类是相同的。这是一个非常简单的电动汽车模型，可以充电和行驶，从而改变其充电状态：

import pandas as pd
import pyomo.environ as pyo


class EV:

    def __init__(self, idx: int, capacity_kwh: float, starting_soc: float,
                 time_steps: pd.DatetimeIndex, driving_profile_km: list[float]):
        self.idx = idx
        self.capacity_kwh = capacity_kwh
        self.soc = starting_soc
        self.driving_profile_km = pd.Series(driving_profile_km, index=time_steps)

    def charge(self, power_kw: float):
        duration_h = 1.0
        self.soc = self.soc + (power_kw * duration_h) / self.capacity_kwh

    def drive(self, distance_km: float):
        consumption_kwh_per_km = 0.15
        self.soc = self.soc - (distance_km * consumption_kwh_per_km) / self.capacity_kwh

两种情况的相同之处在于模型创建的开始以及时间步长和电动车辆的添加：

model = pyo.ConcreteModel()

time_steps = pd.date_range("2023-11-15 10:00", "2023-11-15 13:00", freq="h")
model.time_steps = pyo.Set(initialize=time_steps)

ev1 = EV(1, capacity_kwh=40.0, starting_soc=0.5, time_steps=time_steps,
         driving_profile_km=[20.0, 0.0, 0.0, 20.0])
ev2 = EV(2, capacity_kwh=30.0, starting_soc=0.5, time_steps=time_steps,
         driving_profile_km=[10.0, 5.0, 0.0, 15.0])
ev3 = EV(3, capacity_kwh=20.0, starting_soc=0.5, time_steps=time_steps,
         driving_profile_km=[0.0, 30.0, 25.0, 5.0])
model.electric_vehicles = pyo.Set(initialize=[ev1, ev2, ev3])

现在是上面提到的创建决策变量的时候了。在第一种情况下，我为每辆电动汽车创建一个变量，并使用时间步对其进行索引：

for ev in model.electric_vehicles:
    model.add_component(f"charging_profile_ev{ev.idx}",
                        pyo.Var(model.time_steps, bounds=[0, 11], initialize=0.0))

在第二种情况下，我仅创建一个变量并使用时间步长和电动汽车对其进行索引：

model.add_component(f"charging_profiles",
                    pyo.Var(model.time_steps, model.electric_vehicles,
                            bounds=[0, 11], initialize=0.0))

根据情况，额外需要的表达式、约束和目标的创建也会略有变化。我将简要展示这两种情况，以便 MWE 完整：

案例1：

def expression_soc(m, ts, ev):
    charging_power_kw = m.component(f"charging_profile_ev{ev.idx}")[ts]
    driving_distance_km = ev.driving_profile_km.loc[ts]
    ev.charge(charging_power_kw)
    ev.drive(driving_distance_km)
    return ev.soc
model.add_component(f"soc_expression",
                    pyo.Expression(model.time_steps, model.electric_vehicles,
                                   rule=expression_soc))

def constraint_soc(m, ts, ev):
    soc_value = m.component(f"soc_expression")[ts, ev]
    return 0.5, soc_value, 1.0
model.add_component(f"soc_constraint",
                    pyo.Constraint(model.time_steps, model.electric_vehicles,
                                   rule=constraint_soc))

def objective(m):
    value = 0
    for ev in m.electric_vehicles:
        for ts in m.time_steps:
            value = value + m.component(f"charging_profile_ev{ev.idx}")[ts]
    return value
model.objective = pyo.Objective(rule=objective, sense=pyo.minimize)

results = pyo.SolverFactory("glpk").solve(model)
    model.objective.display()
    model.charging_profile_ev1.display()
    model.charging_profile_ev2.display()
    model.charging_profile_ev3.display()

案例2：

def expression_soc(m, ts, ev):
    charging_power_kw = m.component(f"charging_profiles")[ts, ev]
    driving_distance_km = ev.driving_profile_km.loc[ts]
    ev.charge(charging_power_kw)
    ev.drive(driving_distance_km)
    return ev.soc
model.add_component(f"soc_expression",
                    pyo.Expression(model.time_steps, model.electric_vehicles,
                                   rule=expression_soc))

def constraint_soc(m, ts, ev):
    soc_value = m.component(f"soc_expression")[ts, ev]
    return 0.5, soc_value, 1.0
model.add_component(f"soc_constraint",
                    pyo.Constraint(model.time_steps, model.electric_vehicles,
                                   rule=constraint_soc))

def objective(m):
    value = 0
    for ev in m.electric_vehicles:
        for ts in m.time_steps:
            value = value + m.component(f"charging_profiles")[ts, ev]
    return value
model.objective = pyo.Objective(rule=objective, sense=pyo.minimize)

results = pyo.SolverFactory("glpk").solve(model)
    model.objective.display()
    model.charging_profiles.display()

相应的输出在逻辑上看起来不同，您必须以不同的方式访问配置文件才能在其他应用程序中使用它们。但是是否还有其他差异或原因应首选其中一种或另一种方法？

Answer 1

100％你想要使用第二条路线并索引所有内容，而不是让一堆单例变量滚动。

• 跟踪一切变得更加容易
• 您可以非常轻松地使用 Pyomo 的索引功能来创建“foreach”约束
• 减少拼写错误等风险

就您的模型而言，您有 2 个自然集合：

time

和

vehicle

，我希望大多数变量都由这 2 个集合进行双索引。

首先，然而（这可能是一个很长的帖子来展示这一点......抱歉！）在我看来，你在模型构建中对待

soc

的方式是在玩火。您隐式依赖实例变量

soc

来覆盖所有时间段，但它只是一个变量。它无法同时代表模型所有时期的 soc。你（有点）可以构建你正在做的表达式，但它很脆弱。含义：它正在为模型使用的 SOC 生成正确的表达式，但在执行此操作时会破坏实例变量。一个例子如下所示。请注意，SOC 的表达式是您所期望的，因为它是按顺序构建的 只要您以正确的顺序传递时间步，但是当我再次使用

ev.soc

变量时，它会在“坏”中损坏约束，因为它仍然保留着最后的用法并且它不知道任何关于时间的信息。合理？您可以（如果您小心的话）继续使用您生成的表达式（请参阅下一个约束），因为它在构建时已捕获了正确的表达式。出于类似的原因，您稍后必须再次使用该表达式来计算 SOC。

脆弱的例子：

"""
EV charge model
"""

import pyomo.environ as pyo

class SimpleEV:
    idx = 1 # rolling counter
    def __init__(self):
        self.soc = 10  # initial charge
        self.idx = SimpleEV.idx
        SimpleEV.idx += 1

    def add_juice(self, charge):
        self.soc += charge

    def use_juice(self, charge):
        self.soc -= charge

    # for cleaner printing in pyomo, provide a __repr__
    def __repr__(self):
        return f'ev{self.idx}'

m = pyo.ConcreteModel('ev')

m.T = pyo.Set(initialize=range(3))      # time
m.E = pyo.Set(initialize=[SimpleEV()])  # 1 simple EV

m.charge = pyo.Var(m.E, m.T, domain=pyo.NonNegativeReals, doc='charge added to vehicle e at time t')

m.draws = pyo.Param(m.T, initialize={1:5}, default=0)

def expr_generator(m, e: SimpleEV, t):
    e.add_juice(m.charge[e, t])
    e.use_juice(m.draws[t])
    return e.soc

m.soc_expression = pyo.Expression(m.E, m.T, rule=expr_generator)

# this is NOT going to work using the .soc variable:
def bad_min_charge(m, e: SimpleEV, t):
    return e.soc >= 2
m.bad_min_charge_constraint = pyo.Constraint(m.E, m.T, rule=bad_min_charge)

def ok_min_charge(m, e: SimpleEV, t):
    return m.soc_expression[e, t] >= 2
m.ok_min_charge_constraint  = pyo.Constraint(m.E, m.T, rule=ok_min_charge)

m.pprint()

输出（截断）：

1 Expression Declarations
    soc_expression : Size=3, Index=soc_expression_index
        Key      : Expression
        (ev1, 0) : 10 + charge[ev1,0]
        (ev1, 1) : 10 + charge[ev1,0] + charge[ev1,1] - 5
        (ev1, 2) : 10 + charge[ev1,0] + charge[ev1,1] - 5 + charge[ev1,2]

2 Constraint Declarations
    bad_min_charge_constraint : Size=3, Index=bad_min_charge_constraint_index, Active=True
        Key      : Lower : Body                                                   : Upper : Active
        (ev1, 0) :   2.0 : 10 + charge[ev1,0] + charge[ev1,1] - 5 + charge[ev1,2] :  +Inf :   True
        (ev1, 1) :   2.0 : 10 + charge[ev1,0] + charge[ev1,1] - 5 + charge[ev1,2] :  +Inf :   True
        (ev1, 2) :   2.0 : 10 + charge[ev1,0] + charge[ev1,1] - 5 + charge[ev1,2] :  +Inf :   True
    ok_min_charge_constraint : Size=3, Index=ok_min_charge_constraint_index, Active=True
        Key      : Lower : Body                  : Upper : Active
        (ev1, 0) :   2.0 : soc_expression[ev1,0] :  +Inf :   True
        (ev1, 1) :   2.0 : soc_expression[ev1,1] :  +Inf :   True
        (ev1, 2) :   2.0 : soc_expression[ev1,2] :  +Inf :   True

更好的想法：我会构建类似于下面的模型。一些注意事项：

• 我会将路线与电动车分开，如图所示。我认为它更干净。
• 只要您不乱用

  __hash__

  和

  __eq__

  功能，您就可以将电动汽车直接安全地放入套件中。
• 您应该在代码中添加

  __repr__()

  ，以便输出更加清晰。 （请参阅 python dox 获取指导，我在这里偷工减料）
• 你的语法对于制作组件来说非常笨拙，我认为下面我的语法对于制作集合、约束等来说更干净。
• 下面的构造将所有状态变量（如 SOC）与对象分开，我认为这更干净/更安全。您可以在优化后回填它们，如我在最后所示。

更好的例子：

"""
improved EV charge model
"""

import pyomo.environ as pyo
class DriveProfile:
    def __init__(self, distances):
        self.distances = distances  # {time:km} dictionary

    def get_distance(self, time):
        return self.distances.get(time, None)  # NONE if no travels at this time block
class EV:
    idx = 1  # rolling counter

    def __init__(self, max_soc=100, efficiency=2, route: DriveProfile = None):
        self.initial_soc = 10  # initial charge
        self.max_soc = max_soc  # max charge
        self.discharge_rate = efficiency  # units/km
        self.route = route
        self.idx = EV.idx
        self.soc = None
        EV.idx += 1

    def load_soc(self, soc_dict: dict):
        self.soc = soc_dict

    # for cleaner printing in pyomo, provide a __repr__
    def __repr__(self):
        return f'ev{self.idx}'

# make some fake data...
times = list(range(4))

r1 = DriveProfile({1: 10, 3: 2.4})
r2 = DriveProfile({2: 20})

e1 = EV(route=r1, efficiency=3)
e2 = EV(max_soc=200, route=r2)

# make the model
m = pyo.ConcreteModel('ev')

# SETS
m.T = pyo.Set(initialize=times)     # time
m.E = pyo.Set(initialize=[e1, e2])  # ev's

# VARS
m.soc = pyo.Var(m.E, m.T, domain=pyo.NonNegativeReals, doc='state of charge of E at time T')
m.juice = pyo.Var(m.E, m.T, domain=pyo.NonNegativeReals, doc='juice added to E at time T')

# PARAMS
#   --- all params are internal to the EV objects.  We could capture them here, but not needed. ---

# OBJ:  Minimize the juice delivered
m.obj = pyo.Objective(expr=sum(m.juice[e, t] for e in m.E for t in m.T))


# CONSTRAINTS

# max charge
def max_charge(m, e: EV, t):
    return m.soc[e, t] <= e.max_soc
m.max_charge = pyo.Constraint(m.E, m.T, rule=max_charge)


def soc_balance(m, e: EV, t):
    # find any use (driving)
    use = e.route.get_distance(t)
    discharge = use * e.discharge_rate if use else 0
    if t == m.T.first():
        return m.soc[e, t] == e.initial_soc - discharge + m.juice[e, t]
    return m.soc[e, t] == m.soc[e, m.T.prev(t)] - discharge + m.juice[e, t]
m.soc_balance = pyo.Constraint(m.E, m.T, rule=soc_balance)

m.pprint()

# SOLVE
solver = pyo.SolverFactory('cbc')
soln = solver.solve(m)
print(soln)

# show the juice profile
m.juice.display()

# pass the soc info back to the ev's
print('The state of charge in the vehicles:')
for ev in m.E:
    # we can use items() to get tuples of ( (e, t), soc ) and screen them, then get the .value from the variable...
    res = {time: soc.value for (e, time), soc in m.soc.items() if e == ev}
    ev.load_soc(res)

    print(ev.__repr__(), ev.soc)

输出（截断的第一部分）：

    soc_balance : Size=8, Index=soc_balance_index, Active=True
        Key      : Lower : Body                                                         : Upper : Active
        (ev1, 0) :   0.0 :                             soc[ev1,0] - (10 + juice[ev1,0]) :   0.0 :   True
        (ev1, 1) :   0.0 :                soc[ev1,1] - (soc[ev1,0] - 30 + juice[ev1,1]) :   0.0 :   True
        (ev1, 2) :   0.0 :                     soc[ev1,2] - (soc[ev1,1] + juice[ev1,2]) :   0.0 :   True
        (ev1, 3) :   0.0 : soc[ev1,3] - (soc[ev1,2] - 7.199999999999999 + juice[ev1,3]) :   0.0 :   True
        (ev2, 0) :   0.0 :                             soc[ev2,0] - (10 + juice[ev2,0]) :   0.0 :   True
        (ev2, 1) :   0.0 :                     soc[ev2,1] - (soc[ev2,0] + juice[ev2,1]) :   0.0 :   True
        (ev2, 2) :   0.0 :                soc[ev2,2] - (soc[ev2,1] - 40 + juice[ev2,2]) :   0.0 :   True
        (ev2, 3) :   0.0 :                     soc[ev2,3] - (soc[ev2,2] + juice[ev2,3]) :   0.0 :   True

11 Declarations: T E soc_index soc juice_index juice obj max_charge_index max_charge soc_balance_index soc_balance

Problem: 
- Name: unknown
  Lower bound: 57.2
  Upper bound: 57.2
  Number of objectives: 1
  Number of constraints: 16
  Number of variables: 16
  Number of nonzeros: 3
  Sense: minimize
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.0
  Wallclock time: 0.0
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: None
      Number of created subproblems: None
    Black box: 
      Number of iterations: 2
  Error rc: 0
  Time: 0.007892131805419922
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

juice : juice added to E at time T
    Size=8, Index=juice_index
    Key      : Lower : Value : Upper : Fixed : Stale : Domain
    (ev1, 0) :     0 :  27.2 :  None : False : False : NonNegativeReals
    (ev1, 1) :     0 :   0.0 :  None : False : False : NonNegativeReals
    (ev1, 2) :     0 :   0.0 :  None : False : False : NonNegativeReals
    (ev1, 3) :     0 :   0.0 :  None : False : False : NonNegativeReals
    (ev2, 0) :     0 :  30.0 :  None : False : False : NonNegativeReals
    (ev2, 1) :     0 :   0.0 :  None : False : False : NonNegativeReals
    (ev2, 2) :     0 :   0.0 :  None : False : False : NonNegativeReals
    (ev2, 3) :     0 :   0.0 :  None : False : False : NonNegativeReals
The state of charge in the vehicles:
ev1 {0: 37.2, 1: 7.2, 2: 7.2, 3: 0.0}
ev2 {0: 40.0, 1: 40.0, 2: 0.0, 3: 0.0}