运筹说 第42期 | 算法介绍之运输问题

       本期继续进行运筹学之运输问题算法的讲解，在运输问题中，如何寻找初始可行解以及判断解的最优性是重点的研究问题。通过上期推文的学习，我们知道在求解运输问题初始调运方案时，沃格尔（Vogel）法与西北角法、最小元素法相比，其求解结果往往更接近最优解。在判断一个运输方案是否为最优解时，位势法（对偶变量法）比闭回路法的计算更便捷。

       因此，本期我们将简单回顾一下Vogel法以及位势法的求解步骤，并重点介绍实现这两种方法的Python及MATLAB相关代码，以帮助大家利用工具快速求解运输问题，做到事半功倍。话不多说，我们一起来看看吧！

一、方法介绍

1、寻找初始基可行解—Vogel法

★方法概述

       沃格尔（Vogel）法又称差值法，该方法考虑到，最初按某一最小单位运价优先安排物品调运时，在后续调运过程中却可能不得不采用运费很高的其他供销点，从而使整个运输费用增加。沃格尔法的基本思想是在运价表中分别计算出各行各列的最小单位运价和次小单位运价之差，并称这两个单位运价之差为该销售地或供应地的罚数，然后按照最小单位运价对罚数最大处安排运输。因为若罚数的值很大，说明不按最小运价组织运输就会造成很大的运费损失。

★求解步骤

① 首先计算运输表中每一行和每一列的次小单位运价和最小单位运价之间的差值，分别称为行罚数和列罚数。

② 选取这些罚数中最大者（若存在最大罚数相同的情况，则任选其中一个）所在的行或列的最小单位运价所在的格子，在格子中给其分配尽可能大的运量，划去该行/该列。

③ 在尚未划去的各行或各列中，重复以上步骤，直到最后一个格子也被分配上运量，得到所求运输问题的初始基可行解。

2、判断解的最优性—位势法

★方法概述

       在得到运输问题的初始基可行解后，应对该解做最优性判别。位势法就是用来判断解的最优性的一种方法，其实质是在求解单纯形表中非基变量的检验数。该方法适用于产地和销地较少的运输问题。

★求解步骤

① 在所求得的初始基可行解的表中增加位势列和位势行；

② 由于基变量的检验数等于0，对这组基变量可以构建位势方程组：

       方程组中表示第i行的位势(i=1,2…m)，表示第j列的位势(j=1,2…n)，表示第i行、第j列单元格的运费，这里的和可正、可负也可为零。

       对于特定的调运方案，其可行解中非零变量的个数为(m+n-1)，所以位势方程组含有(m+n-1)个方程。由于运输表中每行和每列都含有基变量，这样构造的方程组含有(m+n)个变量，也就意味着位势的个数为(m+n)。由此可知，位势的个数多于方程的个数，因此在求解时，我们常任意指定某一位势等于一个较小的整数或零，此操作对得到的检验数表无影响。

③ 计算检验数：，得到检验数表。

④ 判断最优性：对于极小化问题来说，如果检验数表中存在<0,说明将变为基变量将使运输费用减少，当前解未达到最优，需要对解进行进一步的调整；否则，当前解达到最优。

3、求解流程总结

       运输问题的求解本质上是一种迭代法，一般的求解步骤为：先计算初始基可行解，再进行解的最优性检验，若得到的解不是最优解，则进行解的改进，并得到新解，再对新解进行最优性检验和改进，直到得到最优解为止。本文选用Vogel法求解初始基可行解、位势法作解的最优性检验，具体的求解流程如下图所示。

二、算法实现

1、例题介绍

       我们依旧沿用上期算法介绍用到的例题，即：某部门有3个生产同类的工厂（产地），生产的产品由4个销售点（销地）出售，各工厂的生产量、各销售点的销售量（假定单位均为t）以及各工厂到各销售点的单位运价（元/t）示于下表中。要求研究产品如何调运才能使总运费最小。

2、平台实现

       此次我们以上述例题为例，借助Python和MATLAB两种平台分别介绍实现寻找初始基可行解的Vogel法以及判断解的最优性的位势法的详细代码、代码调用方法以及最终的运行结果。由于篇幅有限，小编接下来只展示部分代码，小伙伴们可以关注“运筹学”公众号→后台回复“运输问题之Python”及“运输问题之MATLAB”分别获取对应的完整代码。

（1）Python

★准备工作

       在Excel中输入例题中的成本矩阵如下，其中A1：D3为各工厂（产地）到各销售点（销地）的单位运价（元/吨），E1：E3为3个工厂（产地）的产量（吨），A4：D4为4个销售点（销地）的销量（吨）。将上例的运输表保存为EXCEL文件，文件名为“Vogel例题”。

★代码展示

       同求解流程所示步骤相同，程序的运行逻辑也是在利用“Vogel法寻找初始基可行解”后再用“位势法判断解的最优性”。“Vogel法寻找初始基可行解”的部分代码如图一所示，“位势法判断解的最优性”的部分代码如图二所示，在程序中输入文件的保存路径并输出最终结果的代码如图三所示，小伙伴们可以关注“运筹学”公众号→后台回复“运输问题之Python”获取完整代码。

★运行结果

       本期例题用Vogel法求得的初始基可行解是：=12，=4，=8，=2，=14，=8，上述元素为基元素，其他变量的值等于0，为非基元素。该初始基可行解的含义为：A1运往B3的产品数量为12t，A1运往B4的产品数量为4t，A2运往B1的产品数量为8t，A2运往B4的产品数量为2t，A3运往B2的产品数量为14t，A3运往B4的产品数量为8t。当前总成本即运费为244元。

       位势法判断解的最优性”结果如下：当前位势：u1=0，u2=-2，u3=-5，v1=4，v2=10，v3=4，v4=11，并且检验数表中所有的σij≥0，该运输问题已经达到最优解。

★注意事项

       将运输问题的成本矩阵保存到EXCEL后，在程序中输入的文件路径以及文件名应与该文件的实际保存路径以及文件名一致，否则程序会报错。

（2）MATLAB

★代码展示

       MATLAB平台需要分别建立两个函数文件，Vogel函数用于求解运输问题的初始基可行解，部分代码如图一所示，Potential函数为位势法，用于判断解的最优性，部分代码如图二所示，小伙伴们可以关注“运筹学”公众号→后台回复“运输问题之MATLAB”获取完整代码。

★代码调用及运行结果

点击运行“Vogel”文件和“Potential”文件。

       步骤一，在命令行窗口依次“输入运价矩阵→回车→输入产量矩阵→回车→输入销量矩阵→回车→输入Bv = Vogel(A,S,D)→回车”，便可得到初始调运方案结果，结果中NaN为非基元素，正常显示数字的为基元素。可以看出，A1运往B3的产品数量为12t，A1运往B4的产品数量为4t，A2运往B1的产品数量为8t，A2运往B4的产品数量为2t，A3运往B2的产品数量为14t，A3运往B4的产品数量为8t。

       步骤二，在命令行窗口“输入cv = sum(sum(A .* Bv, 'omitnan')) →回车”，可以得到此时的运费为244元。

       步骤三，在命令行窗口“输入SigmaV = Potential(A, Bv) →回车”，可以得到初始调运方案中非基变量的检验数，发现此时检验数均满足大于等于零的条件，由此可以判断步骤一得到的初始基可行解即为最优解。

★注意事项

       在MATLAB命令行窗口输入口令时，输入状态一定要为英文状态，输入的标点符号如果为中文格式的会报错。

• 参考资料：

【Python实现】利用伏格尔 (Vogel) 法寻找初始基可行解

【Python实现】运输问题的表上作业法（一）：利用伏格尔 (Vogel) 法寻找初始基可行解_m0_46927406的博客-CSDN博客_表上作业法伏格尔法

【Python实现】利用位势法判断当前解的最优性

【Python实现】运输问题的表上作业法（二）：利用位势法判断当前解的最优性_m0_46927406的博客-CSDN博客_位势法例题

【Matlab实现】Vogel法求运输问题初始基解

【运筹学】Vogel法求运输问题初始基解的Matlab实现 - 知乎

【Matlab实现】位势法判断运输问题解的最优性

【运筹学】运输问题用位势法计算检验数的Matlab实现 - 知乎

本期的内容就介绍到这里，想要进一步了解运筹学，关注本公众号，快快学起来吧！

作者 |尹萌娟 曹贵玲

责编 | 何洋洋

审核 | 徐小峰

· 知乎｜运筹说 ·

· Bilibili｜运筹说 ·

· CSDN｜运筹说 ·