由于好奇,看了一下代码,然后想研究一下源代码,这个模型最简单直白的理解大体如下:
能够更准确地考虑每种产品的重要性和对运输网络的影响,同时确保运输计划的合理性和经济性。这些因素对于优化解决方案的制定和实施都具有重要的指导意义。
我们来看看代码,对应案例multicommodity_dw.py
#
# This file is part of the Cardinal Optimizer, all rights reserved.
#
import coptpy as cp
from coptpy import COPT
import math
import itertools
# Optimization data for multicommodity problem
ORIG = ['GARY', 'CLEV', 'PITT']
DEST = ['FRA', 'DET', 'LAN', 'WIN', 'STL', 'FRE', 'LAF']
PROD = ['bands', 'coils', 'plate']
supply_list = [400, 800, 200,
700, 1600, 300,
800, 1800, 300]
supply = dict(zip(itertools.product(ORIG, PROD), supply_list))
demand_list = [300, 500, 100,
300, 750, 100,
100, 400, 0,
75, 250, 50,
650, 950, 200,
225, 850, 100,
250, 500, 250]
demand = dict(zip(itertools.product(DEST, PROD), demand_list))
limit_list = [625.0] * len(ORIG) * len(DEST)
limit = dict(zip(itertools.product(ORIG, DEST), limit_list))
cost_list = [30, 39, 41,
10, 14, 15,
8, 11, 12,
10, 14, 16,
11, 16, 17,
71, 82, 86,
6, 8, 8,
22, 27, 29,
7, 9, 9,
10, 12, 13,
7, 9, 9,
21, 26, 28,
82, 95, 99,
13, 17, 18,
19, 24, 26,
11, 14, 14,
12, 17, 17,
10, 13, 13,
25, 28, 31,
83, 99, 104,
15, 20, 20]
cost = dict(zip(itertools.product(ORIG, DEST, PROD), cost_list))
导库加上参数设置,这一步没看懂的同学老师可以看看我昨天发的:
有问题评论区指出吧!
# Optimization parameters for Dantzig-Wolfe decomposition
itercnt = 0
priceconvex = 1.0
price = {(i, j): 0 for i in ORIG for j in DEST}
propcost = []
propshipl = []
vweight = []
现在重头戏来了,就是这个参数设置
这个vweight是运输权重,propship是总运输量约束。
这个方法叫做dw(Dantzig-Wolfe)方法,我们这样理解把,这个方法是将问题分解成master问题和sub问题的。然后迭代优化寻找全局最优解。
迭代优化这个词好好体会一下!
# Create COPT environment
env = cp.Envr()
# Create the MASTER and SUB problem
mMaster = env.createModel()
mSub = env.createModel()
# Disable log information
mMaster.setParam(COPT.Param.Logging, 0)
mSub.setParam(COPT.Param.Logging, 0)
创造环境,然后构建两个子问题的求解环境,这样理解应该没问题。
# Add variable to the MASTER problem
vexcess = mMaster.addVar(name='excess')
#在名为 mMaster 的线性规划模型中添加了一个变量 vexcess,这个变量通常代表MASTER问题中的余量或者剩余量。name='excess' 是变量的名称标识。
# Add constraints to the MASTER problem
cmulti = mMaster.addConstrs((-vexcess <= limit[i, j] for i in ORIG for j in DEST), \
nameprefix='multi')
#向 mMaster 中添加了多个约束 cmulti。每个约束的形式是 -vexcess <= limit[i, j],其中 limit[i, j] 可能是一个参数或者变量,表示限制条件。这些约束适用于所有的原产地 i 和目的地 j 组合。nameprefix='multi' 指定了约束的名称前缀。
# Add initial 'convex' constraint to the MASTER problem
cconvex = mMaster.addConstr(cp.LinExpr(0.0) == 1.0, 'convex')
#向 mMaster 中添加了一个名为 cconvex 的约束,其形式是 cp.LinExpr(0.0) == 1.0,意味着创建了一个恒等于1的线性表达式。这种类型的约束在一些优化算法中可能用于控制问题的凸性或者可行性。
# Add variables to the SUB problem
vtrans = mSub.addVars(ORIG, DEST, PROD, nameprefix='trans')
#在名为 mSub 的线性规划模型中添加了多个变量 vtrans。这些变量代表从每个原产地 i 到每个目的地 j 的每种产品 k 的运输量。nameprefix='trans' 是变量名的前缀。
# Add constraints to the SUB problem
mSub.addConstrs((vtrans.sum(i, '*', k) == supply[i, k] for i in ORIG for k in PROD), \
nameprefix='supply')
mSub.addConstrs((vtrans.sum('*', j, k) == demand[j, k] for j in DEST for k in PROD), \
nameprefix='demand')
#两行代码分别添加了供应约束和需求约束到 mSub 模型中。第一行确保每个原产地 i 的每种产品 k 的总运输量等于该原产地的供应量 supply[i, k]。第二行确保每个目的地 j 的每种产品 k 的总运输量等于该目的地的需求量 demand[j, k]。
# Set objective function for MASTER problem
mMaster.setObjective(vexcess, COPT.MINIMIZE)
#代码设置了 mMaster 模型的优化目标函数。在这个例子中,目标函数是最小化 vexcess,即优化过程将试图尽量减少 vexcess 这个变量的值。
注释我放在代码里面了,可以简单了解一下。
# Dantzig-Wolfe decomposition
print(" *** Dantzig-Wolfe Decomposition *** ")
# Phase I
print('Phase I: ')
while True:
# Display the iteration number
print('Iteration {}: '.format(itercnt))
# Set objective function 'Artificial Reduced Cost' for the SUB problem
mSub.setObjective(cp.quicksum(-price[i, j] * vtrans[i, j, k] for i in ORIG for j in DEST for k in PROD) \
- priceconvex, \
COPT.MINIMIZE)
# Solve the SUB problem
mSub.solve()
做这段算法的过程,迭代优化这样的。
# Check if problem is feasible
if mSub.objval >= -1e-6:
print('No feasible solution...')
break
else:
itercnt += 1
# Get the solution of 'trans' variables in the SUB problem
transval = mSub.getInfo(COPT.Info.Value, vtrans)
# Calculate parameters for the MASTER problem
propship = {(i, j): transval.sum(i, j, '*') for i in ORIG for j in DEST}
propcost.append(transval.prod(cost))
propshipl.append(propship)
# Update constraints 'multi' in the MASTER problem
weightcol = cp.Column()
weightcol.addTerms(cmulti, propship)
weightcol.addTerm(cconvex, 1.0)
# Add variables 'weight' to the MASTER problem
vweight.append(mMaster.addVar(column=weightcol))
# Solve the MASTER problem
mMaster.solve()
# Update 'price'
if mMaster.objval <= 1e-6:
break
else:
price = mMaster.getInfo(COPT.Info.Dual, cmulti)
priceconvex = cconvex.pi
# Phase II
print('Phase II: ')
# Fix variable 'excess'
vexcess_x = vexcess.x
vexcess.lb = vexcess_x
vexcess.ub = vexcess_x
# Set objective function 'Total Cost' for the MASTER problem
mMaster.setObjective(cp.quicksum(propcost[i] * vweight[i] for i in range(itercnt)), \
COPT.MINIMIZE)
# Solve the MASTER problem
mMaster.solve()
# Update 'price' and 'priceconvex'
price = mMaster.getInfo(COPT.Info.Dual, cmulti)
priceconvex = cconvex.pi
# Set the objection function 'Reduced Cost' for the SUB problem
mSub.setObjective(cp.quicksum((cost[i, j, k] - price[i, j]) * vtrans[i, j, k] \
for i in ORIG for j in DEST for k in PROD) - priceconvex, \
COPT.MINIMIZE)
while True:
print('Iteration {}: '.format(itercnt))
# Solve the SUB problem
mSub.solve()
if mSub.objval >= -1e-6:
print('Optimal solution...')
break
else:
itercnt += 1
# Get the solution of 'trans' variables in the SUB problem
transval = mSub.getInfo(COPT.Info.Value, vtrans)
# Calculate parameters for the MASTER problem
propship = {(i, j): transval.sum(i, j, '*') for i in ORIG for j in DEST}
propcost.append(transval.prod(cost))
propshipl.append(propship)
# Update constraints 'multi' in the MASTER problem
weightcol = cp.Column()
weightcol.addTerms(cmulti, propship)
weightcol.addTerm(cconvex, 1.0)
# Add variables 'weight' to the MASTER problem
vweight.append(mMaster.addVar(obj=propcost[itercnt - 1], column=weightcol))
# Solve the MASTER problem
mMaster.solve()
# Update 'price' and 'priceconvex'
price = mMaster.getInfo(COPT.Info.Dual, cmulti)
priceconvex = cconvex.pi
# Set the objection function 'Reduced Cost' for the SUB problem
mSub.setObjective(cp.quicksum((cost[i, j, k] - price[i, j]) * vtrans[i, j, k] \
for i in ORIG for j in DEST for k in PROD) - priceconvex, \
COPT.MINIMIZE)
# Phase III
print('Phase III: ')
optship = {(i, j): 0 for i in ORIG for j in DEST}
for i in ORIG:
for j in DEST:
optship[i, j] = sum(propshipl[k][i, j] * vweight[k].x for k in range(itercnt))
# Set objective function 'Opt Cost' for SUB problem
mSub.setObjective(vtrans.prod(cost), COPT.MINIMIZE)
# Add new constraints to the SUB problem
mSub.addConstrs((vtrans.sum(i, j, '*') == optship[i, j] for i in ORIG for j in DEST))
# Solve the SUB problem
mSub.solve()
print(" *** End Loop *** \n")
# Report solution
print(" *** Summary Report *** ")
print('Objective value: {}'.format(mSub.objval))
print('Variable solution: ')
for i in ORIG:
for j in DEST:
for k in PROD:
if math.fabs(vtrans[i, j, k].x) >= 1e-6:
print(' {0:s} = {1:.6f}'.format(vtrans[i, j, k].name, vtrans[i, j, k].x))
这段代码可以简单理解为分别对两个子问题进行攻击。
Cardinal Optimizer v7.1.3. Build date Apr 29 2024
Copyright Cardinal Operations 2024. All Rights Reserved
Dantzig-Wolfe Decomposition ***
Phase I:
Iteration 0:
Iteration 1:
Iteration 2:
Iteration 3:
Iteration 4:
Phase II:
Iteration 5:
Iteration 6:
Iteration 7:
Iteration 8:
Iteration 9:
Iteration 10:
Iteration 11:
Iteration 12:
Iteration 13:
Iteration 14:
Iteration 15:
Iteration 16:
Iteration 17:
Iteration 18:
Iteration 19:
Iteration 20:
Iteration 21:
Iteration 22:
Iteration 23:
Iteration 24:
Iteration 25:
Optimal solution...
Phase III:
End Loop ***
Summary Report ***
Objective value: 199500.0
Variable solution:
400.000000 =
25.000000 =
200.000000 =
625.000000 =
150.000000 =
225.000000 =
50.000000 =
225.000000 =
300.000000 =
100.000000 =
400.000000 =
173.461464 =
50.000000 =
250.000000 =
300.000000 =
201.538536 =
100.000000 =
225.000000 =
75.000000 =
500.000000 =
50.000000 =
75.000000 =
450.000000 =
100.000000 =
75.000000 =
76.538536 =
625.000000 =
225.000000 =
23.461464 =
250.000000 =
125.000000 =
250.000000 =
求解结果如上。
I, Phase II, Phase III):
Phase I:
0 到 Iteration 4),通过解决MASTER问题和更新相关参数来逐步优化。每次迭代都尝试改进目标函数值,直到满足终止条件或达到一定的迭代次数。
Phase II:
Phase II 继续优化,进行了更多次的迭代(Iteration 5 到 Iteration 25),通过调整变量和约束条件来进一步优化目标函数。最终在Iteration 25 找到了一个最优解。
Phase III:
Phase III 是优化过程的最后阶段。在这个阶段,程序确认找到了一个最优解,并且已经不需要进一步的迭代。因此,打印了 "*** End Loop ***" 表示优化过程结束。
Summary Report:
最后,程序打印了一个汇总报告,展示了优化后的结果:
(Objective value): 199500.0,表示经过优化后的成本或效益。
(Variable solution): 列出了每个关键变量的最优值。例如,trans(GARY, STL, bands) 的运输量为400.0,trans(GARY, FRE, coils) 的运输量为625.0等等。
这段运行结果表明优化过程成功找到了一个满足约束条件的最优解,用于解决特定的多商品运输问题,同时每个阶段和迭代都通过打印消息来展示优化过程的进行情况。
手牵手共建国产科学计算软件和运筹优化软件生态。