集合

声明

可以使用Set和RangeSet类的实例或通过分配集合表达式来声明集合。最简单的集合声明创建一个集合并推迟其成员的创建：

model.A = pyo.Set()

Set 类接受可选参数，例如：

  • dimen = 集合成员的维度

  • doc = 描述集合的字符串

  • filter = 一个布尔函数，在构造过程中用于指示是否应将潜在的新成员分配给集合

  • initialize = 一个包含集合初始成员的可迭代对象，或者返回集合初始成员的可迭代对象的函数。

  • ordered = 一个布尔指示器，表示集合是否有序；默认值为 True

  • validate = 一个布尔函数，用于验证新成员数据

  • within = 用于验证的集合；它是被声明集合的超集。

通常，Pyomo 在构建时会尝试推断集合组件的“维度”（即明显的索引数量）。然而，在某些情况下，Pyomo 可能无法检测到维度（例如，一个未初始化的 Set），或者用户可能希望断言集合的维度。这可以通过 dimen 关键字来实现。例如，要创建一个成员为包含两个项目的元组的集合，可以这样写：

model.B = pyo.Set(dimen=2)

要创建集合 model.A 中所有数字的两倍集合，可以使用

def DoubleA_init(model):
    return (i*2 for i in model.A)
model.C = pyo.Set(initialize=DoubleA_init)

顺便提一下，我们注意到在Python中，总是有很多方法可以完成同样的事情。另外，请注意，如果model.A包含的元素没有定义乘以二的操作，这将会生成一个错误。

initialize 选项可以接受任何 Python 可迭代对象，包括 set、list 或 tuple。这些数据可以从函数返回或直接指定，如下所示：

model.D = pyo.Set(initialize=['red', 'green', 'blue'])

initialize 选项还可以指定一个生成器或函数来定义集合成员。如果是生成器，生成器产生的所有数据将成为初始集合成员：

def X_init(m):
    for i in range(10):
        yield 2*i+1
model.X = pyo.Set(initialize=X_init)

对于初始化函数，Pyomo支持两种签名。第一种是函数返回一个可迭代对象（set、list或tuple），其中包含用于初始化集合的数据：

def Y_init(m):
    return [2*i+1 for i in range(10)]
model.Y = pyo.Set(initialize=Y_init)

在第二个签名中，函数为每个元素调用，将元素编号作为额外参数传递。这个过程会重复，直到函数返回特殊值 Set.End：

def Z_init(model, i):
    if i > 10:
        return pyo.Set.End
    return 2*i+1
model.Z = pyo.Set(initialize=Z_init)

请注意，元素编号从1开始，而不是0：

>>> model.X.pprint()
X : Size=1, Index=None, Ordered=Insertion
    Key  : Dimen : Domain : Size : Members
    None :     1 :    Any :   10 : {1, 3, 5, 7, 9, 11, 13, 15, 17, 19}
>>> model.Y.pprint()
Y : Size=1, Index=None, Ordered=Insertion
    Key  : Dimen : Domain : Size : Members
    None :     1 :    Any :   10 : {1, 3, 5, 7, 9, 11, 13, 15, 17, 19}
>>> model.Z.pprint()
Z : Size=1, Index=None, Ordered=Insertion
    Key  : Dimen : Domain : Size : Members
    None :     1 :    Any :   10 : {3, 5, 7, 9, 11, 13, 15, 17, 19, 21}

有关集合初始化的迭代器的更多信息在 [PyomoBookIII] 书中。

注意

对于抽象模型，在输入文件中指定的数据或通过data参数传递给AbstractModel.create_instance()的数据将覆盖由初始化选项指定的数据。

如果集合作为参数传递给Set而没有关键字，它们将被解释为集合数组的索引。例如，要创建一个由集合model.A的成员索引的集合数组，请使用：

model.E = pyo.Set(model.A)

参数可以组合使用。例如，要创建一个集合数组，按集合 model.A 索引，其中每个集合包含三维成员，请使用：

model.F = pyo.Set(model.A, dimen=3)

initialize 选项可用于创建一个包含一系列数字的集合，但 RangeSet 类为简单序列提供了一种简洁的机制。该类接受起始值、最终值和步长作为其参数。如果 RangeSet 只有一个参数，则该值定义序列中的最终值；起始值和步长默认为1。如果给出两个值，则它们是序列中的第一个和最后一个值，步长默认为1。例如，以下声明创建了一个包含数字1.5、5和8.5的集合：

model.G = pyo.RangeSet(1.5, 10, 3.5)

操作

集合也可以通过使用其他Pyomo集合的集合操作结果来创建。Pyomo支持的集合操作包括并集、交集、差集和对称差集：

model.I = model.A | model.D # union
model.J = model.A & model.D # intersection
model.K = model.A - model.D # difference
model.L = model.A ^ model.D # exclusive-or

例如，叉积运算符是星号（*）。要定义一个新集合 M，它是集合 B 和 C 的叉积，可以使用

model.M = model.B * model.C

这将创建一个虚拟集合，该集合持有对原始集合的引用， 因此对原始集合（B 和 C）的任何更新都会反映 在新集合（M）中。 相比之下，您还可以创建一个具体 集合，它直接存储创建时的笛卡尔积值，并且不会反映原始 集合的后续更改：

model.M_concrete = pyo.Set(initialize=model.B * model.C)

最后，你可以表示集合的成员被限制在两个其他集合的笛卡尔积中，可以使用within关键字：

model.N = pyo.Set(within=model.B * model.C)

预定义虚拟集

为了指定集合、参数和变量的域，Pyomo 提供了以下预定义的虚拟集合：

  • Any = 所有可能的值

  • Reals = 浮点数值

  • PositiveReals = 严格正浮点数值

  • NonPositiveReals = 非正浮点数值

  • NegativeReals = 严格为负的浮点数值

  • NonNegativeReals = 非负浮点数值

  • PercentFraction = 区间 [0,1] 内的浮点数值

  • UnitInterval = PercentFraction 的别名

  • Integers = 整数值

  • PositiveIntegers = 正整数

  • NonPositiveIntegers = 非正整数

  • NegativeIntegers = 负整数值

  • NonNegativeIntegers = 非负整数值

  • Boolean = 布尔值，可以表示为 False/True, 0/1, 'False'/'True' 和 'F'/'T'

  • Binary = 整数 {0, 1}

例如，如果集合 model.O 被声明为在虚拟集合 NegativeIntegers 内，那么尝试添加任何非负整数将导致错误。以下是声明：

model.O = pyo.Set(within=pyo.NegativeIntegers)

稀疏索引集

集合为参数、变量和其他集合提供索引。索引集合的问题对建模者来说很重要，部分原因是出于效率考虑，但主要是因为正确选择索引集合可以产生非常自然的公式，有助于理解和维护。Pyomo利用Python提供了丰富的索引集合创建和使用选项。

如何表示索引的选择通常取决于应用程序和预期的实例数据的性质。为了说明一些选项和问题，我们将考虑涉及网络的问题。在许多网络应用中，声明一组节点是有用的，例如

model.Nodes = pyo.Set()

然后可以参照节点创建一组弧线。

考虑以下最小成本流问题的简单版本：

\begin{array}{lll} \mbox{minimize} & \sum_{a \in \mathcal{A}} c_{a}x_{a} \\ \mbox{subject to:} & S_{n} + \sum_{(i,n) \in \mathcal{A}}x_{(i,n)} & \\ & -D_{n} - \sum_{(n,j) \in \mathcal{A}}x_{(n,j)} & n \in \mathcal{N} \\ & x_{a} \geq 0, & a \in \mathcal{A} \end{array}

其中

  • 集合：节点 \(\equiv \mathcal{N}\)

  • 集合：弧线 \(\equiv \mathcal{A} \subseteq \mathcal{N} \times \mathcal{N}\)

  • 变量：弧上的流量 \((i,j)\) \(\equiv x_{i,j},\; (i,j) \in \mathcal{A}\)

  • 参数: 弧上的流量成本 \((i,j)\) \(\equiv c_{i,j},\; (i,j) \in \mathcal{A}\)

  • 参数：节点latexmath:i \(\equiv D_{i},\; i \in \mathcal{N}\)的需求

  • 参数：节点latexmath:i \(\equiv S_{i},\; i \in \mathcal{N}\)的供应

在最简单的情况下，弧可以只是节点的叉积，这是通过定义实现的

model.Arcs = model.Nodes*model.Nodes

创建一个包含二维成员的集合。对于所有节点始终连接到所有其他节点的应用程序，这可能足够了。然而，当网络不完全密集时，可能会出现一些问题。例如，避免在不存在的弧上流动的负担落在数据文件上，必须为这些弧提供足够高的成本。这样的方案不是很优雅或健壮。

对于许多网络流应用，使用声明弧可能更好

model.Arcs = pyo.Set(dimen=2)

或

model.Arcs = pyo.Set(within=model.Nodes*model.Nodes)

区别在于第一个版本在数据被分配到集合元素时提供错误检查。这将使得稀疏网络的指定更加自然。但这导致需要更改FlowBalance约束，因为在简单示例中编写时，它为每个节点对整个节点集进行求和。一种补救方法是仅对集合model.arcs的成员进行求和，如下所示：

def FlowBalance_rule(m, node):
    return m.Supply[node] \
        + sum(m.Flow[i, node] for i in m.Nodes if (i,node) in m.Arcs) \
        - m.Demand[node] \
        - sum(m.Flow[node, j] for j in m.Nodes if (j,node) in m.Arcs) \
        == 0

除非在稀疏网络中节点数量变得非常大，否则这将没有问题，然后生成此约束的时间可能会成为一个问题（诚然，仅对于非常大的网络，但这样的网络确实存在）。

另一种方法，在许多网络应用中非常有用，是为每个节点设置一个集合，该集合包含指向当前节点的弧的另一端的节点，另一个集合则给出出向弧上的节点。如果这些集合分别称为model.NodesIn和model.NodesOut，那么流量平衡规则可以重写为

def FlowBalance_rule(m, node):
    return m.Supply[node] \
        + sum(m.Flow[i, node] for i in m.NodesIn[node]) \
        - m.Demand[node] \
        - sum(m.Flow[node, j] for j in m.NodesOut[node]) \
        == 0

NodesIn 和 NodesOut 的数据可以添加到输入文件中，这可能是最有效的选择。

对于大多数网络来说，与其从数据文件中读取Arcs、NodesIn和NodesOut，不如只从数据文件中读取Arcs，并通过initialize选项声明model.NodesIn，指定如下创建方式：

def NodesIn_init(m, node):
    for i, j in m.Arcs:
        if j == node:
            yield i
model.NodesIn = pyo.Set(model.Nodes, initialize=NodesIn_init)

与model.NodesOut的定义类似。这段代码为NodesIn创建了一个集合列表，每个节点对应一个节点集合。完整的模型是：

import pyomo.environ as pyo

model = pyo.AbstractModel()

model.Nodes = pyo.Set()
model.Arcs = pyo.Set(dimen=2)

def NodesOut_init(m, node):
    for i, j in m.Arcs:
        if i == node:
            yield j
model.NodesOut = pyo.Set(model.Nodes, initialize=NodesOut_init)

def NodesIn_init(m, node):
    for i, j in m.Arcs:
        if j == node:
            yield i
model.NodesIn = pyo.Set(model.Nodes, initialize=NodesIn_init)

model.Flow = pyo.Var(model.Arcs, domain=pyo.NonNegativeReals)
model.FlowCost = pyo.Param(model.Arcs)

model.Demand = pyo.Param(model.Nodes)
model.Supply = pyo.Param(model.Nodes)

def Obj_rule(m):
    return pyo.summation(m.FlowCost, m.Flow)
model.Obj = pyo.Objective(rule=Obj_rule, sense=pyo.minimize)

def FlowBalance_rule(m, node):
    return m.Supply[node] \
        + sum(m.Flow[i, node] for i in m.NodesIn[node]) \
        - m.Demand[node] \
        - sum(m.Flow[node, j] for j in m.NodesOut[node]) \
        == 0
model.FlowBalance = pyo.Constraint(model.Nodes, rule=FlowBalance_rule)

对于这个模型，一个玩具数据文件（以AMPL的“.dat”格式）将是：

set Nodes := CityA CityB CityC ;

set Arcs :=
CityA CityB
CityA CityC
CityC CityB
;

param : FlowCost :=
CityA CityB 1.4
CityA CityC 2.7
CityC CityB 1.6
 ;

param Demand :=
CityA 0
CityB 1
CityC 1
;

param Supply :=
CityA 2
CityB 0
CityC 0
;

这也可以通过使用BuildAction（更多信息，请参见BuildAction and BuildCheck）来更高效且可能更清晰地完成：

model.NodesOut = pyo.Set(model.Nodes, within=model.Nodes)
model.NodesIn = pyo.Set(model.Nodes, within=model.Nodes)

def Populate_In_and_Out(model):
    # loop over the arcs and record the end points
    for i, j in model.Arcs:
        model.NodesIn[j].add(i)
        model.NodesOut[i].add(j)

model.In_n_Out = pyo.BuildAction(rule=Populate_In_and_Out)

稀疏索引集示例

人们可能希望有一个约束条件成立

\[\forall \; i \in I, k \in K, v \in V_k\]

有许多方法可以完成这个任务，但一个好方法是创建一组由所有model.k, model.V[k]对组成的元组。这可以按如下方式完成：

def kv_init(m):
    return ((k,v) for k in m.K for v in m.V[k])
model.KV = pyo.Set(dimen=2, initialize=kv_init)

我们现在可以创建约束 \(x_{i,k,v} \leq a_{i,k}y_i \;\forall\; i \in I, k \in K, v \in V_k\) 使用：

model.a = pyo.Param(model.I, model.K, default=1)

model.y = pyo.Var(model.I)
model.x = pyo.Var(model.I, model.KV)

def c1_rule(m, i, k, v):
    return m.x[i,k,v] <= m.a[i,k]*m.y[i]
model.c1 = pyo.Constraint(model.I, model.KV, rule=c1_rule)