操作Pyomo模型

本节概述了在使用Pyomo模型时常用的脚本命令。这些命令必须应用于具体模型实例，换句话说，即已实例化的模型。

重复求解

>>> import pyomo.environ as pyo
>>> from pyomo.opt import SolverFactory
>>> model = pyo.ConcreteModel()
>>> model.nVars = pyo.Param(initialize=4)
>>> model.N = pyo.RangeSet(model.nVars)
>>> model.x = pyo.Var(model.N, within=pyo.Binary)
>>> model.obj = pyo.Objective(expr=pyo.summation(model.x))
>>> model.cuts = pyo.ConstraintList()
>>> opt = SolverFactory('glpk')
>>> opt.solve(model) 

>>> # Iterate, adding a cut to exclude the previously found solution
>>> for i in range(5):
...    expr = 0
...    for j in model.x:
...        if pyo.value(model.x[j]) < 0.5:
...            expr += model.x[j]
...        else:
...            expr += (1 - model.x[j])
...    model.cuts.add( expr >= 1 )
...    results = opt.solve(model)
...    print ("\n===== iteration",i)
...    model.display() 

为了说明Pyomo的Python脚本，我们考虑一个例子，该例子在文件iterative1.py中，并使用命令执行

python iterative1.py

注意

这是一个包含Pyomo元素的Python脚本，因此使用python命令执行。可以使用pyomo命令，但当Pyomo完成脚本并尝试将结果发送给求解器时，会出现一些奇怪的消息，这是pyomo命令的功能。

此脚本创建一个模型，解决它，然后添加一个约束以排除刚刚找到的解决方案。这个过程会重复进行，因此脚本会找到并打印多个解决方案。它创建的特定模型只是四个二进制变量的总和。人们不需要计算机来解决这个问题，甚至不需要迭代解决方案。提供此示例只是为了说明脚本编写的一些基本方面。

#  ___________________________________________________________________________
#
#  Pyomo: Python Optimization Modeling Objects
#  Copyright (c) 2008-2024
#  National Technology and Engineering Solutions of Sandia, LLC
#  Under the terms of Contract DE-NA0003525 with National Technology and
#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
#  rights in this software.
#  This software is distributed under the 3-clause BSD License.
#  ___________________________________________________________________________

# iterative1.py
import pyomo.environ as pyo
from pyomo.opt import SolverFactory


# Create a solver
opt = pyo.SolverFactory('glpk')

#
# A simple model with binary variables and
# an empty constraint list.
#
model = pyo.AbstractModel()
model.n = pyo.Param(default=4)
model.x = pyo.Var(pyo.RangeSet(model.n), within=pyo.Binary)


def o_rule(model):
    return pyo.summation(model.x)


model.o = pyo.Objective(rule=o_rule)
model.c = pyo.ConstraintList()

# Create a model instance and optimize
instance = model.create_instance()
results = opt.solve(instance)
instance.display()

# Iterate to eliminate the previously found solution
for i in range(5):
    expr = 0
    for j in instance.x:
        if pyo.value(instance.x[j]) == 0:
            expr += instance.x[j]
        else:
            expr += 1 - instance.x[j]
    instance.c.add(expr >= 1)
    results = opt.solve(instance)
    print("\n===== iteration", i)
    instance.display()

现在让我们分析这个脚本。第一行是一个注释，恰好给出了文件的名称。接下来是两行导入Pyomo符号的代码。pyomo命名空间被导入为pyo。因此，每次使用Pyomo名称时，前面必须加上pyo.。

# iterative1.py
import pyomo.environ as pyo
from pyomo.opt import SolverFactory

通过调用SolverFactory并提供一个指定求解器名称的参数来创建一个执行优化的对象。例如，如果希望使用Gurobi而不是glpk，参数将是'gurobi'：

# Create a solver
opt = pyo.SolverFactory('glpk')

注释后的几行代码创建了一个模型。在我们这里的讨论中，我们将这个模型称为基础模型，因为它将在之后通过添加约束进行扩展。（“基础模型”这个词并不是保留词汇，它们只是为了讨论这个例子而引入的）。基础模型中没有约束，但这只是为了保持简单。基础模型中可能存在约束。尽管它是一个抽象模型，但基础模型通过这些命令完全指定，因为它不需要外部数据：

model = pyo.AbstractModel()
model.n = pyo.Param(default=4)
model.x = pyo.Var(pyo.RangeSet(model.n), within=pyo.Binary)


def o_rule(model):
    return pyo.summation(model.x)


model.o = pyo.Objective(rule=o_rule)

下一行不是基础模型规范的一部分。它创建了一个空的约束列表，脚本将使用该列表来添加约束。

model.c = pyo.ConstraintList()

下一行非注释行创建了实例化的模型，并使用Python变量instance引用实例对象。使用pyomo脚本运行的模型通常不包含这一行，因为模型实例化是由pyomo脚本完成的。在这个例子中，create函数在没有参数的情况下被调用，因为不需要任何参数；然而，在许多脚本中，带有数据命令的文件名是作为参数给出的。

instance = model.create_instance()

下一行调用求解器，并使用Python变量results引用包含结果的对象。

results = opt.solve(instance)

solve函数将结果加载到实例中，因此下一行会写出更新后的值。

instance.display()

下一行非注释行是一个Python迭代命令，它将依次将0到4的整数赋值给Python变量i，尽管该变量在脚本中并未使用。这个循环是导致脚本生成五个更多解决方案的原因：

for i in range(5):

一个表达式在名为 expr 的 Python 变量中构建。Python 变量 j 将被迭代地赋值为变量 x 的所有索引。对于每个索引，变量的值（由刚刚描述的 load 方法加载）将被测试是否为零，并且 expr 中的表达式将相应地增加。尽管 expr 被初始化为 0（一个整数），当它被赋值为涉及 Pyomo 变量对象的表达式时，其类型将变为 Pyomo 表达式：

    expr = 0
    for j in instance.x:
        if pyo.value(instance.x[j]) == 0:
            expr += instance.x[j]
        else:
            expr += 1 - instance.x[j]

在第一次迭代时（当 i 为 0 时），我们知道所有 x 的值都将为 0，因此我们可以预测表达式将是什么样子。我们知道 x 由 1 到 4 的整数索引，因此我们知道 j 将取 1 到 4 的值，并且我们还知道所有索引的 x 值都将为零，因此我们知道 expr 的值将类似于

0 + instance.x[1] + instance.x[2] + instance.x[3] + instance.x[4]

j 的值将被评估，因为它是一个 Python 变量；然而，因为它是一个 Pyomo 变量，instance.x[j] 的值不会被使用，而是变量对象将出现在表达式中。这正是我们在此情况下想要的。当我们在 if 语句中想要使用当前值时，我们使用了 value 函数来获取它。

下一行向名为 c 的约束列表添加了表达式必须大于或等于一的要求：

    instance.c.add(expr >= 1)

证明这排除了最后一种解决方案留给读者作为练习。

外部for循环的最后几行找到解决方案并显示它：

    results = opt.solve(instance)
    print("\n===== iteration", i)
    instance.display()

注意

将求解输出分配给结果对象有些过时。许多脚本只是使用

>>> opt.solve(instance) 

由于结果默认被移动到实例中，结果对象中几乎没有什么有趣的内容。如果出于某种原因，你希望结果保留在结果对象中，不移动到实例中，你可以使用

>>> results = opt.solve(instance, load_solutions=False) 

如果担心求解器没有以最优解终止，这种方法可能有用。例如，

>>> results = opt.solve(instance, load_solutions=False) 
>>> if results.solver.termination_condition == TerminationCondition.optimal: 
...     instance.solutions.load_from(results) 

更改模型或数据并重新求解

上面的iterative1.py示例说明了如何更改模型然后重新求解。在那个示例中，模型通过添加约束来更改，但模型也可以通过更改参数的值来更改。然而，请注意，在这些示例中，我们对具体模型实例进行了更改。这对于AbstractModel用户来说尤为重要，因为这意味着使用instance对象而不是model对象，这使我们能够避免为每次求解创建新的model对象。以下是AbstractModel用户的基本思路：

1. 创建一个AbstractModel（假设它被称为model）

2. 调用 model.create_instance() 来创建一个实例（假设它被称为 instance）

3. 解决 instance

4. 在instance中更改某些内容

5. 再次解决 instance

注意

使用ConcreteModel的用户通常将他们的模型命名为model，这可能会对文档的新手读者造成混淆。基于AbstractModel的示例将引用instance，而使用ConcreteModel的用户通常会使用名称model。

如果 instance 有一个名为 Theta 的参数，该参数被声明为 mutable（即 mutable=True），并且其索引包含 idx，那么可以使用 NewVal 中的值来赋值给它。

>>> instance.Theta[idx] = NewVal

对于一个名为 sigma 的单例参数（即如果它没有被索引），可以使用以下方式进行赋值

>>> instance.sigma = NewVal

注意

如果Param未被声明为可变的，尝试对其进行赋值时将会发生错误。

有关访问Pyomo参数的更多信息，请参阅本文档中关于Param访问的部分访问参数值。请注意，对于具体模型，模型就是实例。

修复变量并重新解决

脚本通常用于固定变量值，而不是更改模型数据。以下示例说明了这一点。

#  ___________________________________________________________________________
#
#  Pyomo: Python Optimization Modeling Objects
#  Copyright (c) 2008-2024
#  National Technology and Engineering Solutions of Sandia, LLC
#  Under the terms of Contract DE-NA0003525 with National Technology and
#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
#  rights in this software.
#  This software is distributed under the 3-clause BSD License.
#  ___________________________________________________________________________

# iterative2.py

import pyomo.environ as pyo
from pyomo.opt import SolverFactory

# Create a solver
opt = pyo.SolverFactory('cplex')

#
# A simple model with binary variables and
# an empty constraint list.
#
model = pyo.AbstractModel()
model.n = pyo.Param(default=4)
model.x = pyo.Var(pyo.RangeSet(model.n), within=pyo.Binary)


def o_rule(model):
    return pyo.summation(model.x)


model.o = pyo.Objective(rule=o_rule)
model.c = pyo.ConstraintList()

# Create a model instance and optimize
instance = model.create_instance()
results = opt.solve(instance)
instance.display()

# "flip" the value of x[2] (it is binary)
# then solve again

if pyo.value(instance.x[2]) == 0:
    instance.x[2].fix(1)
else:
    instance.x[2].fix(0)

results = opt.solve(instance)
instance.display()

在这个例子中，变量是二进制的。模型被求解后，model.x[2]的值被翻转到相反的值，然后再次求解模型。主要关注的行是：

if pyo.value(instance.x[2]) == 0:
    instance.x[2].fix(1)
else:
    instance.x[2].fix(0)

results = opt.solve(instance)

这也可以通过设置上限和下限来实现：

>>> if instance.x[2].value == 0:
...     instance.x[2].setlb(1)
...     instance.x[2].setub(1)
... else:
...     instance.x[2].setlb(0)
...     instance.x[2].setub(0)

请注意，在使用边界时，我们不将fixed设置为True，因为那样会将变量固定在其当前值，然后求解器会忽略边界。

有关访问Pyomo变量的更多信息，请参阅本文档中关于Var访问的部分Accessing Variable Values。

请注意

>>> instance.x.fix(1)

等同于

>>> instance.x.value = 1
>>> instance.x.fixed = True

and

>>> instance.x.fix()

等同于

>>> instance.x.fixed = True

扩展目标函数

可以向ConcreteModel（或实例化的AbstractModel）的目标函数添加项，使用目标函数对象的expr属性。以下是一个简单的示例：

>>> import pyomo.environ as pyo
>>> from pyomo.opt import SolverFactory

>>> model = pyo.ConcreteModel()

>>> model.x = pyo.Var(within=pyo.PositiveReals)
>>> model.y = pyo.Var(within=pyo.PositiveReals)

>>> model.sillybound = pyo.Constraint(expr = model.x + model.y <= 2)

>>> model.obj = pyo.Objective(expr = 20 * model.x)

>>> opt = SolverFactory('glpk') 
>>> opt.solve(model) 

>>> model.pprint() 

>>> print ("------------- extend obj --------------") 
>>> model.obj.expr += 10 * model.y

>>> opt.solve(model) 
>>> model.pprint() 

激活和停用目标

可以声明多个目标，但一次只能激活一个（目前，Pyomo 不支持任何可以接受多个目标的求解器）。如果 model.obj1 和 model.obj2 都已使用 Objective 声明，那么可以确保 model.obj2 被传递给求解器，如下面的简单示例所示：

>>> model = pyo.ConcreteModel()
>>> model.obj1 = pyo.Objective(expr = 0)
>>> model.obj2 = pyo.Objective(expr = 0)

>>> model.obj1.deactivate()
>>> model.obj2.activate()

对于抽象模型，这将在实例化之前完成，否则activate和deactivate调用将针对实例而不是模型。

激活和停用约束

可以使用deactivate()方法暂时禁用约束。 当模型发送到求解器时，不活动的约束不会被包含在内。 可以使用activate()方法重新启用已禁用的约束。

>>> model = pyo.ConcreteModel()
>>> model.v = pyo.Var()
>>> model.con = pyo.Constraint(expr=model.v**2 + model.v >= 3)
>>> model.con.deactivate()
>>> model.con.activate()

索引约束可以整体或按单个索引停用/激活：

>>> model = pyo.ConcreteModel()
>>> model.s = pyo.Set(initialize=[1,2,3])
>>> model.v = pyo.Var(model.s)
>>> def _con(m, s):
...    return m.v[s]**2 + m.v[s] >= 3
>>> model.con = pyo.Constraint(model.s, rule=_con)
>>> model.con.deactivate()   # Deactivate all indices
>>> model.con[1].activate()  # Activate single index