以一道线性规划问题为例,使用Pyomo构建简单的规划模型,调用求解器进行求解。
线性规划、数学建模、Pyomo、cbc求解器
✦
Pyomo简介...
Pyomo(Python Optimization Modeling Objects)是一个用于建模和解决优化问题的Python库。通过声明式建模语言定义规划问题的框架,调用求解器(Gurobi,GLPK,CBC,CPLEX,IPOPT)进行求解。模块化是它的优点,创建和维护起来会更轻松。
✦
过程演示...
from pyomo.environ import *
# 创建一个简单的线性规划模型
model = ConcreteModel()
# 定义变量
model.x = Var(within=NonNegativeReals)
# 定义目标函数和约束
model.obj = Objective(expr=model.x, sense=minimize)
model.con = Constraint(expr=model.x >= 5)
# 选择求解器(例如GLPK)并求解模型
solver = SolverFactory('glpk')
solver.solve(model)
# 打印结果
print("Variable x =", value(model.x))
✦
例题...
在未来的六个月内, 公司计划生产四种元件, 型号分别为 X43-M1, X43-M2, Y54-N1,Y54-N2。这些元件的产量受到产能变化的影响,并且每次产能改变后都需要重新进行控制和调整, 因此会带来不可忽略的费用。因此公司希望最小化这些改变带来的费用,以及生产和库存的成本。
在下表中列出了每种产品每个时期内的需求量,生产和库存成本,初始库存量,以及最后希望保留的库存量。
另外,当产量发生变化时, 需要对设备进行调整。由此带来的费用与当月产量较前一月的改变量(提高或降低)成正比。产量每提高一个产品单位, 则需要多支出 1 欧元;产量每降低一个产品单位, 只需要多支出 0.50 欧元。
目标是最小化由于产量改变引起的费用,以及生产和库存成本,给出总成本及每个月每种产品的生产量。
✦
建模设计一...
集合
月份:Month = [1,2,3,4,5,6]
产品:Products = ['X43-M1','X43-M2','Y54-N1','Y54-N2']
参数
Demands_{i,j}:产品 i 在 j 月的需求量
ProductCost_i:产品 i 的生产成本
StorageCost_i:产品 i 的存储成本
InitialStor_i:产品 i 的期初库存
FinalStor_i:产品 i 的期末库存
定义变量
Production_{i,j}:产品 i 在 j 月的生产量, 非负整数
Storage_{i,j}:产品 i 在 j 月的存储量, 非负整数
Increase_j:j 月比上一月增加的生产量, 非负整数
Decrease_j:j 月比上一月减少的生产量, 非正整数
定义目标函数
成本最小:总生产成本+总存储成本+总调整成本
定义约束
1、生产库存平衡约束
2、生产量调整约束
✦
代码求解...
from pyomo.environ import *
#参数
Demands = [[1500, 3000, 2000, 4000, 2000, 2500],
[1300, 800, 800, 1000, 1100, 900],
[2200, 1500, 2900, 1800, 1200, 2100],
[1400, 1600, 1500, 1000, 1100, 1200]]
ProductCost = [20, 25, 10, 15]
StorageCost = [0.4, 0.5, 0.3, 0.3]
InitialStor = [10, 0, 0, 0]
FinalStor = [50, 10, 10, 10]
#变量维度
model = ConcreteModel()
model.I = Set(initialize = [i for i in range(len(Demands))])
model.J = Set(initialize = [i for i in range(len(Demands[0]))])
model.J1 = Set(initialize = [i for i in range(len(Demands[0]) - 1)])
#变量
model.Production = Var(model.I, model.J, within = NonNegativeIntegers, bounds = (0,10000), initialize=0)
model.Storage = Var(model.I, model.J, within = NonNegativeIntegers, bounds = (0,10000), initialize=0)
model.Increase = Var(model.J1, within = NonNegativeIntegers, bounds = (0,10000), initialize=0)
model.Decrease = Var(model.J1, within = NonPositiveIntegers, bounds = (-10000,0), initialize=0)
#目标函数表达式
model.VarProductCost = Expression(expr=sum(ProductCost[i] * model.Production[i, j] for i in model.I for j in model.J))
model.VarChangeCost = Expression(expr=sum(model.Increase[j] - 0.5 * model.Decrease[j] for j in model.J1))
model.VarStorageCost = Expression(expr=sum(StorageCost[i] * model.Storage[i, j] for i in model.I for j in model.J))
#目标函数
model.cost = Objective(expr=model.VarProductCost + model.VarChangeCost + model.VarStorageCost, sense=minimize)
#生产库存平衡约束
model.StorageCost = ConstraintList()
for i in model.I:
for j in model.J:
if j == 0:
model.StorageCost.add(model.Storage[i, j] == InitialStor[i] + model.Production[i, j] - Demands[i][j])
elif j == len(model.J) - 1:
model.StorageCost.add(model.Storage[i, j] == model.Storage[i, j-1] + model.Production[i, j] - Demands[i][j])
model.StorageCost.add(model.Storage[i, j] == FinalStor[i])
else:
model.StorageCost.add(model.Storage[i, j] == model.Storage[i, j-1] + model.Production[i, j] - Demands[i][j])
#生产量调整约束
def diff_rule(model,j):
return sum(model.Production[i, j + 1] - model.Production[i, j] for i in model.I) == model.Increase[j] + model.Decrease[j]
model.conCapOver = Constraint(model.J1, rule=diff_rule)
#求解
opt = SolverFactory('cbc')
result = opt.solve(model)
print('++++++++++++++++++++++++++++++++++++++')
model.Production.pprint()
print('++++++++++++++++++++++++++++++++++++++')
model.Increase.pprint()
print('++++++++++++++++++++++++++++++++++++++')
model.Decrease.pprint()
print('++++++++++++++++++++++++++++++++++++++')
model.Storage.pprint()
print('Finally Cost: ', model.cost())
result.write()
✦
结果展示...
Production : Size=24, Index=Production_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
(0, 0) : 0 : 1490.0 : 10000 : False : False : NonNegativeIntegers
(0, 1) : 0 : 3000.0 : 10000 : False : False : NonNegativeIntegers
(0, 2) : 0 : 2000.0 : 10000 : False : False : NonNegativeIntegers
(0, 3) : 0 : 4000.0 : 10000 : False : False : NonNegativeIntegers
(0, 4) : 0 : 2000.0 : 10000 : False : False : NonNegativeIntegers
(0, 5) : 0 : 2550.0 : 10000 : False : False : NonNegativeIntegers
(1, 0) : 0 : 1300.0 : 10000 : False : False : NonNegativeIntegers
(1, 1) : 0 : 800.0 : 10000 : False : False : NonNegativeIntegers
(1, 2) : 0 : 800.0 : 10000 : False : False : NonNegativeIntegers
(1, 3) : 0 : 1000.0 : 10000 : False : False : NonNegativeIntegers
(1, 4) : 0 : 1100.0 : 10000 : False : False : NonNegativeIntegers
(1, 5) : 0 : 910.0 : 10000 : False : False : NonNegativeIntegers
(2, 0) : 0 : 2882.0 : 10000 : False : False : NonNegativeIntegers
(2, 1) : 0 : 1672.0 : 10000 : False : False : NonNegativeIntegers
(2, 2) : 0 : 2046.0 : 10000 : False : False : NonNegativeIntegers
(2, 3) : 0 : 1800.0 : 10000 : False : False : NonNegativeIntegers
(2, 4) : 0 : 1890.0 : 10000 : False : False : NonNegativeIntegers
(2, 5) : 0 : 1420.0 : 10000 : False : False : NonNegativeIntegers
(3, 0) : 0 : 1400.0 : 10000 : False : False : NonNegativeIntegers
(3, 1) : 0 : 1600.0 : 10000 : False : False : NonNegativeIntegers
(3, 2) : 0 : 2227.0 : 10000 : False : False : NonNegativeIntegers
(3, 3) : 0 : 273.0 : 10000 : False : False : NonNegativeIntegers
(3, 4) : 0 : 1100.0 : 10000 : False : False : NonNegativeIntegers
(3, 5) : 0 : 1210.0 : 10000 : False : False : NonNegativeIntegers
# ==========================================================
Finally Cost: 684209.4
# ==========================================================
Time: 0.04673957824707031
✦
建模设计二(大M法)...
参数
Demands_{i,j}:产品 i 在 j 月的需求量
ProductCost_i:产品 i 的生产成本
StorageCost_i:产品 i 的存储成本
InitialStor_i:产品 i 的期初库存
FinalStor_i:产品 i 的期末库存
M:极大正数
定义变量
Production_{i,j}:产品 i 在 j 月的生产量, 非负整数
Storage_{i,j}:产品 i 在 j 月的存储量, 非负整数
Change_j:j 月相对上一月调整成本, 非负整数
Sita_j:0-1变量,为1代表 j 月的下一月产量增加,为0减少
定义目标函数
成本最小:总生产成本+总存储成本+总调整成本
定义约束
1、生产库存平衡约束
2、生产量调整约束
✦
代码求解...
from pyomo.environ import *
#参数
Demands = [[1500, 3000, 2000, 4000, 2000, 2500],
[1300, 800, 800, 1000, 1100, 900],
[2200, 1500, 2900, 1800, 1200, 2100],
[1400, 1600, 1500, 1000, 1100, 1200]]
ProductCost = [20, 25, 10, 15]
StorageCost = [0.4, 0.5, 0.3, 0.3]
InitialStor = [10, 0, 0, 0]
FinalStor = [50, 10, 10, 10]
#变量维度
model = ConcreteModel()
model.I = Set(initialize = [i for i in range(len(Demands))])
model.J = Set(initialize = [i for i in range(len(Demands[0]))])
model.J1 = Set(initialize = [i for i in range(len(Demands[0]) - 1)])
#变量
model.Production = Var(model.I, model.J, within = NonNegativeIntegers, bounds = (0,10000), initialize=0)
model.Storage = Var(model.I, model.J, within = NonNegativeIntegers, bounds = (0,10000), initialize=0)
model.Change = Var(model.J1, within = NonNegativeIntegers, bounds = (0,10000), initialize=0)
model.sita = Var(model.J1, within=Binary, initialize=1)
#目标函数表达式
model.VarProductCost = Expression(expr=sum(ProductCost[i] * model.Production[i, j] for i in model.I for j in model.J))
model.VarChangeCost = Expression(expr=sum(model.Change[j] for j in model.J1))
model.VarStorageCost = Expression(expr=sum(StorageCost[i] * model.Storage[i, j] for i in model.I for j in model.J))
#目标函数
model.cost = Objective(expr=model.VarProductCost + model.VarChangeCost + model.VarStorageCost, sense=minimize)
#生产库存平衡约束
model.StorageCost = ConstraintList()
for i in model.I:
for j in model.J:
if j == 0:
model.StorageCost.add(model.Storage[i, j] == InitialStor[i] + model.Production[i, j] - Demands[i][j])
elif j == len(model.J) - 1:
model.StorageCost.add(model.Storage[i, j] == model.Storage[i, j-1] + model.Production[i, j] - Demands[i][j])
model.StorageCost.add(model.Storage[i, j] == FinalStor[i])
else:
model.StorageCost.add(model.Storage[i, j] == model.Storage[i, j-1] + model.Production[i, j] - Demands[i][j])
#生产量调整约束
model.ChangeCost = ConstraintList()
for j in model.J1:
model.ChangeCost.add(sum(model.Production[i, j] - model.Production[i, j + 1] for i in model.I) <= (1 - model.sita[j]) * 9999999)
model.ChangeCost.add(sum(model.Production[i, j] - model.Production[i, j + 1] for i in model.I) >= - model.sita[j] * 9999999)
model.ChangeCost.add(model.Change[j] >= sum(model.Production[i, j + 1] - model.Production[i, j] for i in model.I) - (1 - model.sita[j]) * 9999999)
model.ChangeCost.add(model.Change[j] <= sum(model.Production[i, j + 1] - model.Production[i, j] for i in model.I) + (1 - model.sita[j]) * 9999999)
model.ChangeCost.add(model.Change[j] >= -0.5 * sum(model.Production[i, j + 1] - model.Production[i, j] for i in model.I) - model.sita[j] * 9999999)
model.ChangeCost.add(model.Change[j] <= -0.5 * sum(model.Production[i, j + 1] - model.Production[i, j] for i in model.I) + model.sita[j] * 9999999)
#求解
opt = SolverFactory('cbc')
result = opt.solve(model)
print('++++++++++++++++++++++++++++++++++++++')
model.Production.pprint()
print('++++++++++++++++++++++++++++++++++++++')
model.Change.pprint()
print('++++++++++++++++++++++++++++++++++++++')
model.Storage.pprint()
print('++++++++++++++++++++++++++++++++++++++')
model.sita.pprint()
print('Finally Cost: ', model.cost())
result.write()
✦
结果展示...
Production : Size=24, Index=Production_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
(0, 0) : 0 : 1490.0 : 10000 : False : False : NonNegativeIntegers
(0, 1) : 0 : 3000.0 : 10000 : False : False : NonNegativeIntegers
(0, 2) : 0 : 2000.0 : 10000 : False : False : NonNegativeIntegers
(0, 3) : 0 : 4000.0 : 10000 : False : False : NonNegativeIntegers
(0, 4) : 0 : 2000.0 : 10000 : False : False : NonNegativeIntegers
(0, 5) : 0 : 2550.0 : 10000 : False : False : NonNegativeIntegers
(1, 0) : 0 : 1300.0 : 10000 : False : False : NonNegativeIntegers
(1, 1) : 0 : 800.0 : 10000 : False : False : NonNegativeIntegers
(1, 2) : 0 : 800.0 : 10000 : False : False : NonNegativeIntegers
(1, 3) : 0 : 1000.0 : 10000 : False : False : NonNegativeIntegers
(1, 4) : 0 : 1100.0 : 10000 : False : False : NonNegativeIntegers
(1, 5) : 0 : 910.0 : 10000 : False : False : NonNegativeIntegers
(2, 0) : 0 : 2882.0 : 10000 : False : False : NonNegativeIntegers
(2, 1) : 0 : 1672.0 : 10000 : False : False : NonNegativeIntegers
(2, 2) : 0 : 2773.0 : 10000 : False : False : NonNegativeIntegers
(2, 3) : 0 : 1073.0 : 10000 : False : False : NonNegativeIntegers
(2, 4) : 0 : 1891.0 : 10000 : False : False : NonNegativeIntegers
(2, 5) : 0 : 1419.0 : 10000 : False : False : NonNegativeIntegers
(3, 0) : 0 : 1400.0 : 10000 : False : False : NonNegativeIntegers
(3, 1) : 0 : 1600.0 : 10000 : False : False : NonNegativeIntegers
(3, 2) : 0 : 1500.0 : 10000 : False : False : NonNegativeIntegers
(3, 3) : 0 : 1000.0 : 10000 : False : False : NonNegativeIntegers
(3, 4) : 0 : 1100.0 : 10000 : False : False : NonNegativeIntegers
(3, 5) : 0 : 1210.0 : 10000 : False : False : NonNegativeIntegers
# ==========================================================
Finally Cost: 684210.2
# ==========================================================
Time: 1.277825117111206
参考资料:
无
未完待续