当前位置：首页 > news >正文

Pyomo、PuLP 和 OR-Tools 解决约束优化问题效率对比

news 2025/8/29 8:26:03

假设我们有三种原油类型：

我们的任务是混合 6000 桶 原油，使其满足：

API ≥ 35
硫含量 ≤ 1.0%
任何原油的用量不超过其可用量

在此场景中，我们的目标是最小化总成本。这可以根据需要改为最大化预期利润（成本与市场价格之间的差额），或者其他目标函数，例如最大化随时间变化的稳定性。

1 Pyomo

import pyomo.environ as pyomodel = pyo.ConcreteModel()
crudes = ['A', 'B', 'C']
cost = {'A': 70, 'B': 80, 'C': 65}
api = {'A': 34, 'B': 40, 'C': 30}
sulfur = {'A': 1.2, 'B': 0.5, 'C': 2.0}
avail = {'A': 5000, 'B': 3000, 'C': 4000}
model.crudes = pyo.Set(initialize=crudes)
model.vol = pyo.Var(model.crudes, domain=pyo.NonNegativeReals)
model.cost = pyo.Objective(expr=sum(model.vol[c] * cost[c] for c in crudes), sense=pyo.minimize)
model.total_volume = pyo.Constraint(expr=sum(model.vol[c] for c in crudes) == 6000)
model.sulfur = pyo.Constraint(expr=sum(model.vol[c]*sulfur[c] for c in crudes) <= 6000*1.0)
model.api = pyo.Constraint(expr=sum(model.vol[c]*api[c] for c in crudes) >= 6000*35)
model.avail = pyo.ConstraintList()
for c in crudes:model.avail.add(model.vol[c] <= avail[c])
solver = pyo.SolverFactory('glpk')
solver.solve(model)

优点： 建模语言丰富，易于扩展（非线性、MIP）
缺点： 略显冗长，需要外部求解器（GLPK, CBC）

2 PuLP

from pulp import *model = LpProblem("CrudeBlend", LpMinimize)
vol = LpVariable.dicts("vol", ['A', 'B', 'C'], lowBound=0)
model += lpSum([vol[i] * cost[i] for i in crudes])
model += lpSum([vol[i] for i in crudes]) == 6000
model += lpSum([vol[i] * sulfur[i] for i in crudes]) <= 6000 * 1.0
model += lpSum([vol[i] * api[i] for i in crudes]) >= 6000 * 35
for i in crudes:model += vol[i] <= avail[i]
model.solve()

优点： 语法简洁，易于学习，内置 CBC 求解器
缺点： 对于非线性或大型模型功能较弱

3 OR-Tools

from ortools.linear_solver import pywraplpsolver = pywraplp.Solver.CreateSolver('GLOP')
vol = {i: solver.NumVar(0, avail[i], i) for i in crudes}
solver.Add(solver.Sum([vol[i] for i in crudes]) == 6000)
solver.Add(solver.Sum([vol[i]*sulfur[i] for i in crudes]) <= 6000 * 1.0)
solver.Add(solver.Sum([vol[i]*api[i] for i in crudes]) >= 6000 * 35)
solver.Minimize(solver.Sum([vol[i] * cost[i] for i in crudes]))
status = solver.Solve()

优点： 速度极快，适用于 Google 规模系统的生产环境
缺点： 语法灵活性较低，对非线性或符号模型的支持较弱

4 推荐

对于小型、易读的线性规划问题，使用 PuLP。
当建模更复杂、多阶段或非线性系统时，使用 Pyomo。
当需要速度、性能或路径规划时，使用 OR-Tools。

这是一个并排测试。问题很简单，但对于这个简单的问题，运行时差异之大令人惊讶。

import pyomo.environ as pyo
from pulp import *
from ortools.linear_solver import pywraplp
import time
import pandas as pdcrudes = ['A', 'B', 'C']
cost = {'A': 70, 'B': 80, 'C': 65}
api = {'A': 34, 'B': 40, 'C': 30}
sulfur = {'A': 1.2, 'B': 0.5, 'C': 2.0}
avail = {'A': 5000, 'B': 3000, 'C': 4000}
target_vol = 6000
sulfur_max = 1.0
api_min = 35results = {}start = time.time()
model_pyo = pyo.ConcreteModel()
model_pyo.crudes = pyo.Set(initialize=crudes)
model_pyo.vol = pyo.Var(model_pyo.crudes, domain=pyo.NonNegativeReals)
model_pyo.cost = pyo.Objective(expr=sum(model_pyo.vol[c] * cost[c] for c in crudes), sense=pyo.minimize)
model_pyo.total_volume = pyo.Constraint(expr=sum(model_pyo.vol[c] for c in crudes) == target_vol)
model_pyo.sulfur = pyo.Constraint(expr=sum(model_pyo.vol[c]*sulfur[c] for c in crudes) <= target_vol*sulfur_max)
model_pyo.api = pyo.Constraint(expr=sum(model_pyo.vol[c]*api[c] for c in crudes) >= target_vol*api_min)
model_pyo.avail = pyo.ConstraintList()
for c in crudes:model_pyo.avail.add(model_pyo.vol[c] <= avail[c])
solver = pyo.SolverFactory('glpk')
solver.solve(model_pyo)
end = time.time()
results['Pyomo'] = {'runtime_sec': round(end - start, 6),'volumes': {c: round(pyo.value(model_pyo.vol[c]), 2) for c in crudes},'cost': round(sum(pyo.value(model_pyo.vol[c]) * cost[c] for c in crudes), 2)
}start = time.time()
model_pulp = LpProblem("CrudeBlend", LpMinimize)
vol_pulp = LpVariable.dicts("vol", crudes, lowBound=0)
model_pulp += lpSum([vol_pulp[i] * cost[i] for i in crudes])
model_pulp += lpSum([vol_pulp[i] for i in crudes]) == target_vol
model_pulp += lpSum([vol_pulp[i] * sulfur[i] for i in crudes]) <= target_vol * sulfur_max
model_pulp += lpSum([vol_pulp[i] * api[i] for i in crudes]) >= target_vol * api_min
for i in crudes:model_pulp += vol_pulp[i] <= avail[i]
model_pulp.solve()
end = time.time()
results['PuLP'] = {'runtime_sec': round(end - start, 6),'volumes': {c: round(vol_pulp[c].varValue, 2) for c in crudes},'cost': round(value(model_pulp.objective), 2)
}start = time.time()
solver = pywraplp.Solver.CreateSolver('GLOP')
vol_ort = {i: solver.NumVar(0, avail[i], i) for i in crudes}
solver.Add(solver.Sum([vol_ort[i] for i in crudes]) == target_vol)
solver.Add(solver.Sum([vol_ort[i]*sulfur[i] for i in crudes]) <= target_vol * sulfur_max)
solver.Add(solver.Sum([vol_ort[i]*api[i] for i in crudes]) >= target_vol * api_min)
solver.Minimize(solver.Sum([vol_ort[i] * cost[i] for i in crudes]))
status = solver.Solve()
end = time.time()
results['OR-Tools'] = {'runtime_sec': round(end - start, 6),'volumes': {c: round(vol_ort[c].solution_value(), 2) for c in crudes},'cost': round(sum(vol_ort[c].solution_value() * cost[c] for c in crudes), 2)
}df = pd.DataFrame(results).T