杉数求解器-多商品问题加强版本

文摘   2024-07-26 11:27   湖北  

由于好奇,看了一下代码,然后想研究一下源代码,这个模型最简单直白的理解大体如下:

能够更准确地考虑每种产品的重要性和对运输网络的影响,同时确保运输计划的合理性和经济性。这些因素对于优化解决方案的制定和实施都具有重要的指导意义。

我们来看看代码,对应案例multicommodity_dw.py

## This file is part of the Cardinal Optimizer, all rights reserved.#
import coptpy as cpfrom coptpy import COPT
import mathimport itertools
# Optimization data for multicommodity problemORIG = ['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))

导库加上参数设置,这一步没看懂的同学老师可以看看我昨天发的:

杉数求解器-多商品问题的求解(数模玩家必看)

杉数求解器—尝试求解2024年亚太赛中文赛C题

有问题评论区指出吧!

# Optimization parameters for Dantzig-Wolfe decompositionitercnt     = 0priceconvex = 1.0price       = {(i, j): 0 for i in ORIG for j in DEST}propcost    = []propshipl   = []
vweight = []

现在重头戏来了,就是这个参数设置

这个vweight是运输权重,propship是总运输量约束。

这个方法叫做dw(Dantzig-Wolfe)方法,我们这样理解把,这个方法是将问题分解成master问题和sub问题的。然后迭代优化寻找全局最优解。

迭代优化这个词好好体会一下!

# Create COPT environmentenv = cp.Envr()
# Create the MASTER and SUB problemmMaster = env.createModel()mSub = env.createModel()
# Disable log informationmMaster.setParam(COPT.Param.Logging, 0)mSub.setParam(COPT.Param.Logging, 0)

创造环境,然后构建两个子问题的求解环境,这样理解应该没问题。

# Add variable to the MASTER problemvexcess = mMaster.addVar(name='excess')#在名为 mMaster 的线性规划模型中添加了一个变量 vexcess,这个变量通常代表MASTER问题中的余量或者剩余量。name='excess' 是变量的名称标识。

# Add constraints to the MASTER problemcmulti = 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 problemcconvex = mMaster.addConstr(cp.LinExpr(0.0) == 1.0, 'convex')#向 mMaster 中添加了一个名为 cconvex 的约束,其形式是 cp.LinExpr(0.0) == 1.0,意味着创建了一个恒等于1的线性表达式。这种类型的约束在一些优化算法中可能用于控制问题的凸性或者可行性。
# Add variables to the SUB problemvtrans = mSub.addVars(ORIG, DEST, PROD, nameprefix='trans')#在名为 mSub 的线性规划模型中添加了多个变量 vtrans。这些变量代表从每个原产地 i 到每个目的地 j 的每种产品 k 的运输量。nameprefix='trans' 是变量名的前缀。
# Add constraints to the SUB problemmSub.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 problemmMaster.setObjective(vexcess, COPT.MINIMIZE)#代码设置了 mMaster 模型的优化目标函数。在这个例子中,目标函数是最小化 vexcess,即优化过程将试图尽量减少 vexcess 这个变量的值。

注释我放在代码里面了,可以简单了解一下。

# Dantzig-Wolfe decompositionprint("               *** Dantzig-Wolfe Decomposition ***               ")
# Phase Iprint('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 IIprint('Phase II: ')
# Fix variable 'excess'vexcess_x = vexcess.xvexcess.lb = vexcess_xvexcess.ub = vexcess_x
# Set objective function 'Total Cost' for the MASTER problemmMaster.setObjective(cp.quicksum(propcost[i] * vweight[i] for i in range(itercnt)), \ COPT.MINIMIZE)
# Solve the MASTER problemmMaster.solve()
# Update 'price' and 'priceconvex'price = mMaster.getInfo(COPT.Info.Dual, cmulti)priceconvex = cconvex.pi
# Set the objection function 'Reduced Cost' for the SUB problemmSub.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 IIIprint('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 problemmSub.setObjective(vtrans.prod(cost), COPT.MINIMIZE)
# Add new constraints to the SUB problemmSub.addConstrs((vtrans.sum(i, j, '*') == optship[i, j] for i in ORIG for j in DEST))
# Solve the SUB problemmSub.solve()print(" *** End Loop *** \n")
# Report solutionprint(" *** 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 2024Copyright 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.0Variable solution: trans(GARY,STL,bands) = 400.000000 trans(GARY,STL,coils) = 25.000000 trans(GARY,STL,plate) = 200.000000 trans(GARY,FRE,coils) = 625.000000 trans(GARY,LAF,coils) = 150.000000 trans(CLEV,FRA,bands) = 225.000000 trans(CLEV,FRA,plate) = 50.000000 trans(CLEV,DET,bands) = 225.000000 trans(CLEV,DET,coils) = 300.000000 trans(CLEV,DET,plate) = 100.000000 trans(CLEV,LAN,coils) = 400.000000 trans(CLEV,WIN,coils) = 173.461464 trans(CLEV,WIN,plate) = 50.000000 trans(CLEV,STL,bands) = 250.000000 trans(CLEV,STL,coils) = 300.000000 trans(CLEV,FRE,coils) = 201.538536 trans(CLEV,FRE,plate) = 100.000000 trans(CLEV,LAF,coils) = 225.000000 trans(PITT,FRA,bands) = 75.000000 trans(PITT,FRA,coils) = 500.000000 trans(PITT,FRA,plate) = 50.000000 trans(PITT,DET,bands) = 75.000000 trans(PITT,DET,coils) = 450.000000 trans(PITT,LAN,bands) = 100.000000 trans(PITT,WIN,bands) = 75.000000 trans(PITT,WIN,coils) = 76.538536 trans(PITT,STL,coils) = 625.000000 trans(PITT,FRE,bands) = 225.000000 trans(PITT,FRE,coils) = 23.461464 trans(PITT,LAF,bands) = 250.000000 trans(PITT,LAF,coils) = 125.000000 trans(PITT,LAF,plate) = 250.000000

求解结果如上。

这段运行结果展示了一个使用了Dantzig-Wolfe分解方法的优化过程,分为三个阶段(Phase I, Phase II, Phase III):
Phase I:在这个阶段中,代码进行了多次迭代(Iteration 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等等。这段运行结果表明优化过程成功找到了一个满足约束条件的最优解,用于解决特定的多商品运输问题,同时每个阶段和迭代都通过打印消息来展示优化过程的进行情况。
害,到这里终于快速把这段代码看完了。相关的案例我放一个:杉数求解器-可行性松弛
这个案例比较难,源代码可以申请试用杉数求解器然后拿到example。如果后面接触了相关内容,看看能不能再讲的容易理解一点。
说直接一点,比原来的问题整了一个小活就是加了一个权重,然后利用dw,将问题分解分别攻击,最后获得答案。

手牵手共建国产科学计算软件和运筹优化软件生态。

师苑数模
发布数模协会培训推文,讲解数模算法。赛题讲解及比赛通知。学校竞赛结果及学校竞赛成绩发布等文章。
 最新文章