本期繼續(xù)進(jìn)行運(yùn)籌學(xué)之運(yùn)輸問題算法的講解，在運(yùn)輸問題中，如何尋找初始可行解以及判斷解的最優(yōu)性是重點(diǎn)的研究問題。通過上期推文的學(xué)習(xí)，我們知道在求解運(yùn)輸問題初始調(diào)運(yùn)方案時，沃格爾（Vogel）法與西北角法、最小元素法相比，其求解結(jié)果往往更接近最優(yōu)解。在判斷一個運(yùn)輸方案是否為最優(yōu)解時，位勢法（對偶變量法）比閉回路法的計(jì)算更便捷。

因此，本期我們將簡單回顧一下Vogel法以及位勢法的求解步驟，并重點(diǎn)介紹實(shí)現(xiàn)這兩種方法的Python及MATLAB相關(guān)代碼，以幫助大家利用工具快速求解運(yùn)輸問題，做到事半功倍。話不多說，我們一起來看看吧！

一、方法介紹

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

★方法概述

沃格爾（Vogel）法又稱差值法，該方法考慮到，最初按某一最小單位運(yùn)價優(yōu)先安排物品調(diào)運(yùn)時，在后續(xù)調(diào)運(yùn)過程中卻可能不得不采用運(yùn)費(fèi)很高的其他供銷點(diǎn)，從而使整個運(yùn)輸費(fèi)用增加。沃格爾法的基本思想是在運(yùn)價表中分別計(jì)算出各行各列的最小單位運(yùn)價和次小單位運(yùn)價之差，并稱這兩個單位運(yùn)價之差為該銷售地或供應(yīng)地的罰數(shù)，然后按照最小單位運(yùn)價對罰數(shù)最大處安排運(yùn)輸。因?yàn)槿袅P數(shù)的值很大，說明不按最小運(yùn)價組織運(yùn)輸就會造成很大的運(yùn)費(fèi)損失。

★求解步驟

① 首先計(jì)算運(yùn)輸表中每一行和每一列的次小單位運(yùn)價和最小單位運(yùn)價之間的差值，分別稱為行罰數(shù)和列罰數(shù)。

② 選取這些罰數(shù)中最大者（若存在最大罰數(shù)相同的情況，則任選其中一個）所在的行或列的最小單位運(yùn)價所在的格子，在格子中給其分配盡可能大的運(yùn)量，劃去該行/該列。

③ 在尚未劃去的各行或各列中，重復(fù)以上步驟，直到最后一個格子也被分配上運(yùn)量，得到所求運(yùn)輸問題的初始基可行解。

2、判斷解的最優(yōu)性—位勢法

★方法概述

在得到運(yùn)輸問題的初始基可行解后，應(yīng)對該解做最優(yōu)性判別。位勢法就是用來判斷解的最優(yōu)性的一種方法，其實(shí)質(zhì)是在求解單純形表中非基變量的檢驗(yàn)數(shù)。該方法適用于產(chǎn)地和銷地較少的運(yùn)輸問題。

★求解步驟

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

② 由于基變量的檢驗(yàn)數(shù)等于0，對這組基變量可以構(gòu)建位勢方程組：

方程組中表示第i行的位勢(i=1,2…m)，表示第j列的位勢(j=1,2…n)，表示第i行、第j列單元格的運(yùn)費(fèi)，這里的和可正、可負(fù)也可為零。

對于特定的調(diào)運(yùn)方案，其可行解中非零變量的個數(shù)為(m+n-1)，所以位勢方程組含有(m+n-1)個方程。由于運(yùn)輸表中每行和每列都含有基變量，這樣構(gòu)造的方程組含有(m+n)個變量，也就意味著位勢的個數(shù)為(m+n)。由此可知，位勢的個數(shù)多于方程的個數(shù)，因此在求解時，我們常任意指定某一位勢等于一個較小的整數(shù)或零，此操作對得到的檢驗(yàn)數(shù)表無影響。

③ 計(jì)算檢驗(yàn)數(shù)，得到檢驗(yàn)數(shù)表。

④ 判斷最優(yōu)性：對于極小化問題來說，如果檢驗(yàn)數(shù)表中存在<0,說明將變?yōu)榛兞繉⑹惯\(yùn)輸費(fèi)用減少，當(dāng)前解未達(dá)到最優(yōu)，需要對解進(jìn)行進(jìn)一步的調(diào)整；否則，當(dāng)前解達(dá)到最優(yōu)。

3、求解流程總結(jié)

運(yùn)輸問題的求解本質(zhì)上是一種迭代法，一般的求解步驟為：先計(jì)算初始基可行解，再進(jìn)行解的最優(yōu)性檢驗(yàn)，若得到的解不是最優(yōu)解，則進(jìn)行解的改進(jìn)，并得到新解，再對新解進(jìn)行最優(yōu)性檢驗(yàn)和改進(jìn)，直到得到最優(yōu)解為止。本文選用Vogel法求解初始基可行解、位勢法作解的最優(yōu)性檢驗(yàn)，具體的求解流程如下圖所示。

二、算法實(shí)現(xiàn)

1、例題介紹

我們依舊沿用上期算法介紹用到的例題，即：某部門有3個生產(chǎn)同類的工廠（產(chǎn)地），生產(chǎn)的產(chǎn)品由4個銷售點(diǎn)（銷地）出售，各工廠的生產(chǎn)量、各銷售點(diǎn)的銷售量（假定單位均為t）以及各工廠到各銷售點(diǎn)的單位運(yùn)價（元/t）示于下表中。要求研究產(chǎn)品如何調(diào)運(yùn)才能使總運(yùn)費(fèi)最小。

2、平臺實(shí)現(xiàn)

此次我們以上述例題為例，借助Python和MATLAB兩種平臺分別介紹實(shí)現(xiàn)尋找初始基可行解的Vogel法以及判斷解的最優(yōu)性的位勢法的詳細(xì)代碼、代碼調(diào)用方法以及最終的運(yùn)行結(jié)果。由于篇幅有限，小編接下來只展示部分代碼，小伙伴們可以關(guān)注微信“運(yùn)籌學(xué)”公眾號→后臺回復(fù)“運(yùn)輸問題之Python”及“運(yùn)輸問題之MATLAB”分別獲取對應(yīng)的完整代碼。

（1）Python

★準(zhǔn)備工作

在Excel中輸入例題中的成本矩陣如下，其中A1：D3為各工廠（產(chǎn)地）到各銷售點(diǎn)（銷地）的單位運(yùn)價（元/噸），E1：E3為3個工廠（產(chǎn)地）的產(chǎn)量（噸），A4：D4為4個銷售點(diǎn)（銷地）的銷量（噸）。將上例的運(yùn)輸表保存為EXCEL文件，文件名為“Vogel例題”。

★代碼展示

同求解流程所示步驟相同，程序的運(yùn)行邏輯也是在利用“Vogel法尋找初始基可行解”后再用“位勢法判斷解的最優(yōu)性”?！癡ogel法尋找初始基可行解”的部分代碼如圖一所示，“位勢法判斷解的最優(yōu)性”的部分代碼如圖二所示，在程序中輸入文件的保存路徑并輸出最終結(jié)果的代碼如圖三所示，小伙伴們可以關(guān)注“運(yùn)籌學(xué)”公眾號→后臺回復(fù)“運(yùn)輸問題之Python”獲取完整代碼。

★運(yùn)行結(jié)果

本期例題用Vogel法求得的初始基可行解是：=12，?=4，=8，=2，=14，=8，上述元素為基元素，其他變量的值等于0，為非基元素。該初始基可行解的含義為：A1運(yùn)往B3的產(chǎn)品數(shù)量為12t，A1運(yùn)往B4的產(chǎn)品數(shù)量為4t，A2運(yùn)往B1的產(chǎn)品數(shù)量為8t，A2運(yùn)往B4的產(chǎn)品數(shù)量為2t，A3運(yùn)往B2的產(chǎn)品數(shù)量為14t，A3運(yùn)往B4的產(chǎn)品數(shù)量為8t。當(dāng)前總成本即運(yùn)費(fèi)為244元。

“位勢法判斷解的最優(yōu)性”結(jié)果如下：當(dāng)前位勢：=0，=-2，=-5，=4，=10，=4，=11，并且檢驗(yàn)數(shù)表中所有的≥0，該運(yùn)輸問題已經(jīng)達(dá)到最優(yōu)解。

★注意事項(xiàng)

將運(yùn)輸問題的成本矩陣保存到EXCEL后，在程序中輸入的文件路徑以及文件名應(yīng)與該文件的實(shí)際保存路徑以及文件名一致，否則程序會報(bào)錯。

（2）MATLAB

★代碼展示

MATLAB平臺需要分別建立兩個函數(shù)文件，Vogel函數(shù)用于求解運(yùn)輸問題的初始基可行解，部分代碼如圖一所示，Potential函數(shù)為位勢法，用于判斷解的最優(yōu)性，部分代碼如圖二所示，小伙伴們可以關(guān)注“運(yùn)籌學(xué)”公眾號→后臺回復(fù)“運(yùn)輸問題之MATLAB”獲取完整代碼。

★代碼調(diào)用及運(yùn)行結(jié)果

點(diǎn)擊運(yùn)行“Vogel”文件和“Potential”文件。

步驟一，在命令行窗口依次“輸入運(yùn)價矩陣→回車→輸入產(chǎn)量矩陣→回車→輸入銷量矩陣→回車→輸入Bv = Vogel(A,S,D)→回車”，便可得到初始調(diào)運(yùn)方案結(jié)果，結(jié)果中NaN為非基元素，正常顯示數(shù)字的為基元素?？梢钥闯?，A1運(yùn)往B3的產(chǎn)品數(shù)量為12t，A1運(yùn)往B4的產(chǎn)品數(shù)量為4t，A2運(yùn)往B1的產(chǎn)品數(shù)量為8t，A2運(yùn)往B4的產(chǎn)品數(shù)量為2t，A3運(yùn)往B2的產(chǎn)品數(shù)量為14t，A3運(yùn)往B4的產(chǎn)品數(shù)量為8t。

步驟二，在命令行窗口“輸入cv = sum(sum(A .* Bv, 'omitnan')) →回車”，可以得到此時的運(yùn)費(fèi)為244元。

步驟三，在命令行窗口“輸入SigmaV = Potential(A, Bv) →回車”，可以得到初始調(diào)運(yùn)方案中非基變量的檢驗(yàn)數(shù)，發(fā)現(xiàn)此時檢驗(yàn)數(shù)均滿足大于等于零的條件，由此可以判斷步驟一得到的初始基可行解即為最優(yōu)解。

★注意事項(xiàng)

在MATLAB命令行窗口輸入口令時，輸入狀態(tài)一定要為英文狀態(tài)，輸入的標(biāo)點(diǎn)符號如果為中文格式的會報(bào)錯。

☆?? 參考資料：

【Python實(shí)現(xiàn)】利用伏格爾 (Vogel) 法尋找初始基可行解

https://blog.csdn.net/m0_46927406/article/details/109903843?spm=1001.2014.3001.5501

【Python實(shí)現(xiàn)】利用位勢法判斷當(dāng)前解的最優(yōu)性

https://blog.csdn.net/m0_46927406/article/details/110070154?spm=1001.2014.3001.5501

【Matlab實(shí)現(xiàn)】Vogel法求運(yùn)輸問題初始基解

https://zhuanlan.zhihu.com/p/62017151

【Matlab實(shí)現(xiàn)】位勢法判斷運(yùn)輸問題解的最優(yōu)性

https://zhuanlan.zhihu.com/p/62034448

本期的內(nèi)容就介紹到這里，想要進(jìn)一步了解運(yùn)籌學(xué)，關(guān)注本公眾號，快快學(xué)起來吧！

作者 |尹萌娟 曹貴玲

責(zé)編 | 何洋洋

審核 | 徐小峰