【非线性求解工具】pyomo进阶使用手册

Pyomo模型的组件

集合(Sets)

如何建立集合

集合可以使用Set和RangeSet类来声明，也可以通过指定集合表达式来声明。最简单的集合声明创建了一个集合，并推迟了其成员的创建。

model.A = pyo.Set()

可设置的参数如下：

• dimen = 集合元素的维度
• doc = 描述这个集合的文字，类型为str
• filter = 一个布尔函数，在构建过程中用于指示一个潜在的新成员是否应该被分配到集合中。
• initialize = 一个包含Set初始成员的迭代器，或者返回·Set·初始成员的迭代器的函数。
• ordered = 一个布尔指标，表明该集合是有序的;默认是True
• validate = 一个验证新成员数据的布尔函数
• within = 用于验证的集合；它是被声明的集合的超集合。

下面对上面的参数进行说明：
一般来说，Pyomo试图在构造Set组件时推断其 “维度”（即明显的索引数）。然而，在某些情况下，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)

initialize可以接受任何 Python 可迭代数据，包括set,list,tuple。这个数据可以从一个函数中返回，也可以直接指定，如在

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

初始化选项也可以指定一个生成器或一个函数来指定集合的成员。在生成器的情况下，生成器产生的所有数据将成为初始集合成员。

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

对于初始化函数，Pyomo支持两种方式。在第一种方式中，函数返回一个包含数据的迭代器（set、list或tuple），用它来初始化Set。

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)

如果集合作为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集合来存储集合操作的结果来创建。包括集合的并(union)，交(intersection),差(difference)和对称差(symmetric difference)

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 = 所有值
• Real = 浮点值（实数）
• PositiveReals = 正浮点数值
• NonPositiveReals = 非正浮点值
• NegativeReals = 负浮点数
• NonNegativeReals = 非负浮点数
• PercentFraction = 区间[0,1]的浮点数
• UnitInterval = 百分数
• Integers = 整数
• PositiveIntegers = 正整数
• NonPositiveIntegers = 非正整数
• NegativeIntegers = 负整数
• NonNegativeIntegers = 非负整数
• Boolean = False/True, 0/1, ’False’/’True’ and ’F’/’T’
• Binary = {0,1}
  例如，如果集合model.O被声明为在虚拟集合NegativeIntegers中，那么试图添加任何非负整数的东西都会导致错误。下面是这个声明。

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

稀疏的索引集

索引集问题对建模者很重要，部分原因是出于效率的考虑，但主要是因为正确的索引集选择可以导致非常自然的公式，有利于理解和维护。Pyomo利用Python为索引集的创建和使用提供了丰富的选项集合。
比如可能想要建立如下的约束：
有很多方法来实现这一点，一个好的方法是创建一个由所有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}\le a_{i,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)

Parameters

"参数 "这个词在很多场合都会用到。当讨论Pyomo模型时，我们用这个词来指代必须提供的数据，以便找到决策变量的最佳（或良好）赋值。参数被声明为Param类的实例，它接受的参数与Set类有些类似。例如，下面的代码片段声明了集合model.A和model.B，然后是一个由model.A和model.B索引的参数model.P

model.A = pyo.RangeSet(1,3)
model.B = pyo.Set()
model.P = pyo.Param(model.A, model.B)

除了作为索引的集合外，Param还接受以下选项：

• default =默认值
• doc = 一个描述参数的字符串。
• initialize = 一个函数（或Python对象），返回用于初始化参数值的数据。
• mutable = 布尔值，表示在Param被初始化后是否允许改变Param的值。
• validate = 一个回调函数，接收模型、建议值和建议值的索引；如果值是有效的，返回True。返回False将产生一个异常。
• within = 用于验证的集合；它规定了有效参数值的域。

这些选项的执行方式与Set的执行方式相同。
例如，给定model.A的值为{1，2，3}，那么有很多方法可以创建一个参数，代表一个主对角线上有9，16，25，其他地方有0的正方形矩阵，这里有两种方法可以做到。首先使用一个Python对象进行初始化。

v={}
v[1,1] = 9
v[2,2] = 16
v[3,3] = 25
model.S1 = pyo.Param(model.A, model.A, initialize=v, default=0)

现在使用一个初始化函数，为每个索引元组自动调用一次
(记住，我们假设model.A包含{1, 2, 3})

def s_init(model, i, j):
    if i == j:
        return i*i
    else:
        return 0.0
model.S2 = pyo.Param(model.A, model.A, initialize=s_init)

在这个例子中，索引集包含了整数，但是索引集不需要是数字。使用字符串是很常见的。
参数值可以通过一个验证函数来检查。在下面的例子中，参数T（由model.A索引）的每个值被检查为大于3.14159。如果提供的值小于这个值，模型的即时运算将被终止，并发出错误信息。编写验证函数时，如果数据有效，应返回True，否则返回False。

t_data = {1: 10, 2: 3, 3: 20}

def t_validate(model, v, i):
    return v > 3.14159

model.T = pyo.Param(model.A, validate=t_validate, initialize=t_data)

这个例子将产生以下错误，表明为T[2]提供的值没有通过验证:

Traceback (most recent call last):
  ...
ValueError: Invalid parameter value: T[2] = '3', value type=<class 'int'>.
    Value failed parameter validation rule

Variables

可选的指令包括：

• bounds=一个函数(或Python对象)，它给出了变量的(下限,上限)边界对。
• domain=一个集合，它是该变量可以取值的超集合。
• initialize=一个函数（或Python对象），给出变量的起始值；这对非线性模型特别重要。
• within=domain
  下面的代码片段说明了这些选项的某些方面，它声明了一个名为model.LumberJack的singleton（即未索引的）变量，它将采用0到6之间的实值，并且初始化为1.5。

model.LumberJack = Var(within=NonNegativeReals, bounds=(0,6), initialize=1.5)

有时可以用Python的赋值语句来代替初始化选项，如：

model.LumberJack = 1.5

对于有索引的变量，界限和初始值通常由一个规则(Python函数)来指定，该规则本身可以引用参数或其他数据。这些规则的正式参数以模型开始，然后是索引。这在下面的代码片断中得到说明，它利用了被声明为lb和ub的Python字典，这些字典被一个函数用来提供边界。

model.A = Set(initialize=['Scones', 'Tea'])
lb = {'Scones':2, 'Tea':4}
ub = {'Scones':5, 'Tea':7}
def fb(model, i):
   return (lb[i], ub[i])
model.PriceToCharge = Var(model.A, domain=PositiveIntegers, bounds=fb)

Objectives

目标是一个变量的函数，它返回一个优化问题试图最大化或最小化的值。Pyomo中的Objective函数可以声明一个目标。下面是这样一个非常简单的函数版本，假设model.x之前已经被声明为一个Var:

def ObjRule(model):
    return 2*model.x[1] + 3*model.x[2]
model.obj1 = pyo.Objective(rule=ObjRule)

Constraints

大多数约束是使用等式或不等式来指定的。例如，如果变量model.x有 "butter "和 "scones’ "这两个索引，那么这个约束条件限制了这些索引的总和必须正好是3。

def teaOKrule(model):
    return(model.x['butter'] + model.x['scones'] == 3)
model.TeaConst = Constraint(rule=teaOKrule)

除了等式和不等式，还可以用三元组的形式(lb, expr, ub) ，意思是lb<=expr<=ub（lb和ub可以是None）。比如说以下两个是同一个意思：

model.x = Var()

def aRule(model):
   return model.x >= 2
model.Boundx = Constraint(rule=aRule)

def bRule(model):
   return (2, model.x, None)
model.boundx = Constraint(rule=bRule)