pyomo

Pyomo - 数学优化建模

Pyomo 允许你使用自然的数学语法（集合、参数、变量、约束）来定义优化问题。它将模型与求解器分离，因此无需更改代码，同一个模型就能由不同的引擎求解。

首先：验证前置条件

pip install pyomo
# Also install a solver (e.g., GLPK for linear/integer problems)
# Conda: conda install -c conda-forge glpk ipopt

适用场景

• 战略规划：长期资源分配或投资规划。
• 流程工程：优化化工厂或炼油厂的运营（非线性）。
• 能源系统：电网调度和机组组合问题。
• 供应链优化：多周期、多商品流问题。
• 非线性规划（NLP）：约束或目标涉及平滑曲线时（例如 x²、log(x)）。
• 随机规划：对优化中的不确定性进行建模。
• 自定义求解器集成：需要使用 IPOPT、SCIP 或 Baron 等特定求解器时。

参考文档

官方文档：http://www.pyomo.org/
GitHub：https://github.com/Pyomo/pyomo
搜索模式：pyo.ConcreteModel、pyo.Constraint、pyo.Objective、pyo.SolverFactory

核心原则

Concrete 模型与 Abstract 模型

• ConcreteModel：在构建模型时定义数据（在 Python/数据科学中最常见）。
• AbstractModel：先定义结构，之后再加载数据（大型工业模型的标准做法）。

组件

• Var：求解器需要求出的未知量。
• Set/Param：定义问题实例的数据。
• Objective：需要最小化或最大化的函数。
• Constraint：变量必须遵循的规则。

求解器

Pyomo 没有自己的求解器。它需要系统中安装外部求解器（例如用于 LP/MIP 的 glpk 或用于 NLP 的 ipopt）。

快速参考

安装

pip install pyomo
# Also install a solver (e.g., GLPK for linear/integer problems)
# Conda: conda install -c conda-forge glpk ipopt

标准导入

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

基本模式 - 简单线性模型（Concrete）

import pyomo.environ as pyo

# 1. Create Model
model = pyo.ConcreteModel()

# 2. Define Variables
model.x = pyo.Var(within=pyo.NonNegativeReals)
model.y = pyo.Var(within=pyo.NonNegativeReals)

# 3. Define Objective (Minimize x + 2*y)
model.obj = pyo.Objective(expr=model.x + 2*model.y, sense=pyo.minimize)

# 4. Define Constraints
model.con1 = pyo.Constraint(expr=3*model.x + 4*model.y >= 12)
model.con2 = pyo.Constraint(expr=2*model.x + 5*model.y >= 10)

# 5. Solve
solver = pyo.SolverFactory('glpk')
results = solver.solve(model)

# 6. Access Results
print(f"x = {pyo.value(model.x)}, y = {pyo.value(model.y)}")

关键规则

✅ 应该做

• 缩放问题 - 与所有优化问题一样，如果变量和约束的数量级相近（例如 0.1 到 100），Pyomo 的效果最佳。
• 使用 pyo.value() - 求解后，始终用 pyo.value(model.x) 包装变量，以获取其数值结果。
• 检查求解器状态 - 使用 results.solver.termination_condition 验证求解器是否得到了最优解。
• 使用 pyo.Set 和 pyo.Param - 对于复杂模型，使用这些组件组织数据，而不是直接使用 Python 列表/字典。
• 指定求解器选项 - 使用 solver.options['tm_limit'] = 60（或根据求解器设置类似选项）来管理执行时间。
• 使用 model.dual - 如果求解器支持，请启用 dual 后缀以获取影子价格（这对经济分析很重要）。

❌ 不应该做

• 不要使用 Python 数学函数 - 在表达式中使用 pyo.log()、pyo.exp()、pyo.sqrt()，而不是 math.log() 或 np.log()。
• 不要在约束中使用 if 语句 - 约束必须是代数表达式。对于条件逻辑，请使用“Big-M”表示法或逻辑重构。
• 避免“难处理”的非线性 - abs(x) 或 max(x, y) 等函数并非光滑函数，可能导致许多求解器无法正常工作。请使用重构形式。
• 不要忘记截距/常数项 - x + y = 10 这样的约束没有问题，但要确保单位一致。

反模式（绝对不要）

import pyomo.environ as pyo
import numpy as np

# ❌ BAD: Using NumPy/Math functions in constraints
# model.con = pyo.Constraint(expr=np.sin(model.x) <= 0.5)

# ✅ GOOD: Use Pyomo-compatible functions
model.con = pyo.Constraint(expr=pyo.sin(model.x) <= 0.5)

# ❌ BAD: Using Python IF for conditional constraints
# if model.x > 10:
#     model.con = pyo.Constraint(expr=model.y <= 5)

# ✅ GOOD: Using Big-M notation (for binary variable z)
# y <= 5 + M * (1 - z)
# x >= 10 - M * (1 - z)

# ❌ BAD: Printing a variable directly
# print(model.x) # Returns a reference object, not a number!

# ✅ GOOD: Use value()
print(pyo.value(model.x))

使用集合和索引建模

构建大型模型

model = pyo.ConcreteModel()

# Data
products = ['A', 'B', 'C']
profit = {'A': 10, 'B': 20, 'C': 15}
limit = 100

# Components
model.P = pyo.Set(initialize=products)
model.x = pyo.Var(model.P, within=pyo.NonNegativeReals)

# Indexed Objective
def obj_rule(model):
    return sum(profit[p] * model.x[p] for p in model.P)
model.obj = pyo.Objective(rule=obj_rule, sense=pyo.maximize)

# Indexed Constraint
def limit_rule(model):
    return sum(model.x[p] for p in model.P) <= limit
model.con = pyo.Constraint(rule=limit_rule)

非线性规划（NLP）

使用 IPOPT 求解器

# Solving: Minimize (x-2)^2 + (y-2)^2
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=0) # Initialization is CRUCIAL for NLP
model.y = pyo.Var(initialize=0)

model.obj = pyo.Objective(expr=(model.x - 2)**2 + (model.y - 2)**2)

# Constraint: x^2 + y <= 1
model.con = pyo.Constraint(expr=model.x**2 + model.y <= 1)

# Solve with IPOPT
solver = pyo.SolverFactory('ipopt')
solver.solve(model)

混合整数线性规划（MILP）

离散决策

# Binary variable: 1 if we open a warehouse, 0 otherwise
model.use_warehouse = pyo.Var(within=pyo.Binary)

# Integer variable: Number of trucks to buy
model.num_trucks = pyo.Var(within=pyo.NonNegativeIntegers)

# Conditional logic: If warehouse is not used, trucks must be 0
# trucks <= Capacity * use_warehouse
model.cap_con = pyo.Constraint(expr=model.num_trucks <= 100 * model.use_warehouse)

实用工作流

1. 饮食问题（经典线性规划）

def solve_diet(foods, nutrients, costs, requirements):
    model = pyo.ConcreteModel()
    model.F = pyo.Set(initialize=foods)
    model.N = pyo.Set(initialize=nutrients)
    
    model.x = pyo.Var(model.F, within=pyo.NonNegativeReals)
    
    # Minimize cost
    model.obj = pyo.Objective(expr=sum(costs[f] * model.x[f] for f in model.F))
    
    # Meet nutrient requirements
    def nutrient_rule(model, n):
        return sum(nutrients[f][n] * model.x[f] for f in model.F) >= requirements[n]
    model.con = pyo.Constraint(model.N, rule=nutrient_rule)
    
    pyo.SolverFactory('glpk').solve(model)
    return {f: pyo.value(model.x[f]) for f in model.F}

2. 混合问题（化工/流程工程）

# Balancing component fractions in a mixture
# Note: Often becomes non-linear (NLP) if both flow and fraction are variables
def blend_optimization(inputs, target_purity):
    model = pyo.ConcreteModel()
    # ... model setup ...
    # con: sum(flow[i] * purity[i]) / sum(flow[i]) == target_purity
    # becomes: sum(flow[i] * purity[i]) == target_purity * sum(flow[i]) (Linearized)

性能优化

求解器选择

• GLPK/CBC：适用于免费的开源 LP/MIP。
• Gurobi/CPLEX：工业标准（求解大型 MIP 时速度极快）。
• IPOPT：最适合平滑非线性（NLP）。

热启动

对于迭代优化，请使用上一次的解作为起点。

# For NLP solvers like IPOPT
model.x.set_value(prev_x_value)
solver.solve(model)

常见陷阱与解决方案

终止条件

务必检查求解器停止的原因。

from pyomo.opt import TerminationCondition

results = solver.solve(model)
if results.solver.termination_condition == TerminationCondition.optimal:
    print("Success")
elif results.solver.termination_condition == TerminationCondition.infeasible:
    print("Check your constraints!")

NLP 中的非凸性

如果 NLP 模型有多个局部最小值，IPOPT 可能会陷入其中。

# ✅ Solution: 
# 1. Provide multiple different initial guesses (multistart).
# 2. Use a global solver like BARON or SCIP.

使用变量进行索引

model.x[model.y] 中的 y 是 Var，这种写法不合法。

# ✅ Solution: Use model.AddElement or binary variable reformulations.

最佳实践

1. 始终对问题进行缩放——变量和约束的数量级应相近
2. 求解后使用 pyo.value() 提取变量的数值结果
3. 检查求解器的终止条件，以验证最优性或诊断问题
4. 使用 pyo.Set 和 pyo.Param 来组织复杂模型中的数据
5. 为 NLP 问题初始化变量——良好的起点至关重要
6. 在表达式中使用兼容 Pyomo 的函数（pyo.log()、pyo.exp() 等），而不是 NumPy/math 函数
7. 使用 Big-M 表示法或二进制变量重构条件逻辑
8. 避免使用 abs() 或 max() 等非光滑函数——改用重构形式
9. 在需要进行经济分析时，启用对偶后缀以获取影子价格
10. 根据问题规模合理设置求解器的时间限制和选项

Pyomo 是将高层次数学抽象转化为已求解的商业和科学问题的终极工具。它能够弥合代数建模与高性能求解器之间的差距，因此成为现代规范性分析的基础。