1、蔬菜運輸問題2012年8月22日摘要本文運用floyd算法求出各蔬菜采購點到每個菜市場的最短運輸距離，然后用lingo軟件計算蔬菜調運費用及預期短缺損失最小的調運方案，緊接著根據(jù)題目要求對算法加以修改得出每個市場短缺率都小于20%的最優(yōu)調運方案，并求出了最佳的供應改進方案。關鍵詞最短路問題 floyd算法 運輸問題一、問題重述光明市是一個人口不到15萬人的小城市。根據(jù)該市的蔬菜種植情況，分別在花市（A），城鄉(xiāng)路口（B）和下塘街（C）設三個收購點，再由各收購點分送到全市的8個菜市場，該市道路情況，各路段距離（單位：100m）及各收購點，菜市場的具體位置見圖1，按常年情況，A,B,C三個收購點每天

2、收購量分別為200，170和160（單位：100 kg）,各菜市場的每天需求量及發(fā)生供應短缺時帶來的損失（元/100kg）見表1.設從收購點至各菜市場蔬菜調運費為1元/(100kg.100m). 7 5 4 8 3 7 A 7 6 B 6 8 5 5 4 7 11 7 4 7 5 6 6 3 5 8 6 6 10 C 10 5 11 圖1表1菜市場每天需求（100 kg）短缺損失（元/100kg）7510608805701010010558905808(a) 為該市設計一個從收購點至個菜市場的定點供應方案，使用于蔬菜調運及預期的短缺損失為最小；(b) 若規(guī)定各菜市場短缺量一律不超過需求量的20

3、%，重新設計定點供應方案(c) 為滿足城市居民的蔬菜供應，光明市的領導規(guī)劃增加蔬菜種植面積，試問增產(chǎn)的蔬菜每天應分別向A,B,C三個采購點供應多少最經(jīng)濟合理。二、問題分析求總的運費最低，可以先求出各采購點到菜市場的最小運費，由于單位重量運費和距離成正比，題目所給的圖1里包含了部分菜市場、中轉點以及收購點之間的距離，（a）題可以用求最短路的方法求出各采購點到菜市場的最短路徑，乘上單位重量單位距離費用就是單位重量各運輸線路的費用，然后用線性方法即可解得相應的最小調運費用及預期短缺損失。第二問規(guī)定各菜市場短缺量一律不超過需求量的20%，只需要在上題基礎上加上新的限制條件，即可得出新的調運方案。第三問

4、可以在第二問的基礎上用靈敏度分析進行求解，也可以建立新的線性問題進行求解。三、模型假設1、各個菜市場、中轉點以及收購點都可以作為中轉點；2、各個菜市場、中轉點以及收購點都可以的最大容納量為610噸；3、假設只考慮運輸費用和短缺費用，不考慮裝卸等其它費用；4、假設運輸?shù)氖卟寺吠局袥]有損耗；5、忽略從種菜場地到收購點的運輸費用。四、符號說明A收購點分送到全市的8個菜市場的供應量分別為a1,b1,c1,d1,e1,f1,g1,h1,B收購點分送到全市的8個菜市場的供應量分別為a2,b2,c2,d2,e2,f2,g2,h2,C收購點分送到全市的8個菜市場的供應量分別為a3,b3,c3,d3,e3,f3

5、,g3,h3,8個菜市場的短缺損失量分別為a,b,c,d,e,f,g,h(單位均為100kg)。五、模型的建立和求解按照問題的分析，首先就要求解各采購點到菜市場的最短距離，在圖論里面關于最短路問題比較常用的是Dijkstra算法，Dijkstra算法提供了從網(wǎng)絡圖中某一點到其他點的最短距離。主要特點是以起始點為中心向外層層擴展，直到擴展到終點為止。但由于它遍歷計算的節(jié)點很多，所以效率較低，實際問題中往往要求網(wǎng)絡中任意兩點之間的最短路距離。如果仍然采用Dijkstra算法對各點分別計算，就顯得很麻煩。所以就可以使用網(wǎng)絡各點之間的矩陣計算法，即Floyd算法。Floyd算法的基本是：從任意節(jié)點i到

6、任意節(jié)點j的最短路徑不外乎2種可能，1是直接從i到j，2是從i經(jīng)過若干個節(jié)點k到j。i到j的最短距離不外乎存在經(jīng)過i和j之間的k和不經(jīng)過k兩種可能，所以可以令k=1，2，3，.，n(n是城市的數(shù)目)，在檢查d(i,j)和d(i,k)+d(k,j)的值；在此d(i,k)和d(k,j)分別是目前為止所知道的i到k和k到j的最短距離。因此d(i,k)+d(k,j)就是i到j經(jīng)過k的最短距離。所以，若有d(i,j)>d(i,k)+d(k,j)，就表示從i出發(fā)經(jīng)過k再到j的距離要比原來的i到j距離短，自然把i到j的d(i,j)重寫為d(i,k)+d(k,j)，每當一個k查完了，d(i,j)就是目前

7、的i到j的最短距離。重復這一過程，最后當查完所有的k時，d(i,j)里面存放的就是i到j之間的最短距離了。對于上面的問題，只要能列出它的距離的鄰接矩陣，如表2所示。便能用floyd法進行計算了。表2 各點的鄰接矩陣圖首先對圖1中的4個純中轉點進行編號，為（9）（12），并標注在圖1中，A、B、C的編號分別為1、14、15，其余點在矩陣中的編號如表1第一行所示。由于數(shù)據(jù)比較復雜，用普通的計算很困難，所以我們可以用MATLAB軟件來編程求解。算法思路：采用標號作業(yè)法,每次迭代產(chǎn)生一個永久標號, 從而生長一顆以V0為根的最短路樹,在這顆樹上每個頂點和根節(jié)點之間的路徑皆為最短路徑。當然此問題也需要表2

8、的各點鄰接矩陣。接下來可以得到Dijkstra法的MATLAB語言。M程序如下：function min,path=dijkstra(w,start,terminal)n=size(w,1); label(start)=0; f(start)=start;for i=1:n if i=start label(i)=inf;end, ends(1)=start; u=start;while length(s)<n for i=1:n ins=0; for j=1:length(s) if i=s(j) ins=1; end, end if ins=0 v=i; if label(v)>

9、;(label(u)+w(u,v) label(v)=(label(u)+w(u,v); f(v)=u; end, end, end v1=0; k=inf; for i=1:n ins=0; for j=1:length(s) if i=s(j) ins=1; end, end if ins=0 v=i; if k>label(v) k=label(v); v1=v; end, end, end s(length(s)+1)=v1; u=v1;endmin=label(terminal); path(1)=terminal;i=1; while path(i)=start path(i

10、+1)=f(path(i); i=i+1 ;end path(i)=start;L=length(path);path=path(L:-1:1);圖2 Dijkstra算法的MATLAB運行圖 如圖2，紅色框內(nèi)的是Dijkstra算法的腳本程序，藍色框內(nèi)的試運行結果，根據(jù)程序顯示s點到j點的最短路就是4百米。在程序中顯示所走的路線為path=1 2，對應點即為直接從A點到達點。但是由于運用dijkstra編程每次都要修改起始點和終點，所以在大部分情況下效率都不高。由于dijkstra算法較為復雜，所以可用Floyd算法用于求最短路徑。Floyd算法思路本質很簡單，即用三個for循環(huán)的嵌套，代碼

11、思路如下：for ( int i = 0; i < 節(jié)點個數(shù); +i )for ( int j = 0; j < 節(jié)點個數(shù); +j )for ( int k = 0; k < 節(jié)點個數(shù); +k )if ( Disik + Diskj < Disij )/ 找到更短路徑Disij = Disik + Diskj;但是這里我們要注意循環(huán)的嵌套順序，如果把檢查所有節(jié)點X放在最內(nèi)層，那么結果將是不正確的，因為這樣便過早的把i到j的最短路徑確定下來了，所以正確的應該是這樣的： for ( k = 0; k < 節(jié)點個數(shù); +k )for ( i = 0; i < 節(jié)點

12、個數(shù); +i )for ( j = 0; j < 節(jié)點個數(shù); +j )if ( Disik + Diskj < Disij ) 找到更短路徑Disij = Disik + Diskj;.那么接下來的問題就是，我們?nèi)绾握页鲎疃搪窂?這里需要借助一個輔助數(shù)組Path，它是這樣使用的：Path(AB)的值如果為P，則表示A節(jié)點到B節(jié)點的最短路徑是A->.->P->B。這樣一來，假設我們要找A->B的最短路徑，那么就依次查找，假設Path(AB)的值為P，那么接著查找Path(AP)，假設Path(AP)的值為L，那么接著查找Path(AL)，假設Path(AL)的

13、值為A，則查找結束，最短路徑為A->L->P->B。那么，如何填充Path的值呢？很簡單，當我們發(fā)現(xiàn)Dis(AX) + Dis(XB) < Dis(AB)成立時，就要把最短路徑改為A->.->X->.->B，而此時，Path(XB)的值是已知的，所以，Path(AB) = Path(XB)。Floyd算法直接在圖的帶權鄰接矩陣中用插入頂點的方法依次遞推地構造出n個矩陣D(1), D(2), , D(n), D(n)是圖的距離矩陣, 同時引入一個后繼點矩陣記錄兩點間的最短路徑。d(i,j) : i到j的距離; path(i,j): i到j的路徑上i

14、的后繼點;輸入帶權鄰接矩陣a(i,j)。賦初值，對所有i,j, d(i,j)¬a(i,j) , path(i,j)¬j,k=l。然后更新d(i,j) , path(i,j)對所有i,j,若d(i,k)+d(k,j)<d(i,j)則d(i,j)=d(i,k)+d(k,j) , path(i,j)path(i,k) , k =k+1。重復上述步驟直到k=n+1。接下來就是代入具體的數(shù)據(jù)了，這里的圖以及鄰接矩陣依舊是修改后的圖1及表2。Floyd算法的MATLAB代碼如下：M程序function D,path,min1,path1=floyd(a,start,termina

15、l)D=a;n=size(D,1);path=zeros(n,n);for i=1:n for j=1:n if D(i,j)=inf path(i,j)=j;end, end, endfor k=1:n for i=1:n for j=1:n if D(i,k)+D(k,j)<D(i,j) D(i,j)=D(i,k)+D(k,j); path(i,j)=path(i,k);end, end, end,endif nargin=3 min1=D(start,terminal); m(1)=start; i=1; path1= ; while path(m(i),terminal)=ter

16、minal k=i+1; m(k)=path(m(i),terminal); i=i+1; end m(i+1)=terminal; path1=m;end 腳本程序的代碼如下：a=0488infinf6infinf7inf4infinfinf407infinfinf5infinfinfinfinfinfinfinf870infinfinfinfinfinf3infinfinf7inf8infinf0inf5infinfinf5inf467infinfinfinfinf0infinfinfinfinfinfinf5infinfinfinfinf5inf0infinfinfinf673inf66

17、5infinfinfinf0infinfinf75infinfinfinfinfinfinfinfinfinf011inf10infinfinf5infinfinfinfinfinfinf110infinfinf6inf107inf35infinfinfinfinf0infinfinf6infinfinfinfinfinf6710infinf0infinfinf84infinf4inf75infinfinfinf0infinfinfinfinfinf653infinf6infinfinf011infinfinf77infinfinfinfinf6infinf110infinfinfinfinf

18、inf6inf510inf8infinfinf0;D, path=floyd(a)運行結果如下：圖3 返回矩陣D圖4 返回矩陣path D(i,j)表示i到j的最短距離; path(i,j)表示i和j之間的最短路徑上頂點i的后繼點。根據(jù)圖3、圖4的結果可以很快的知道各點到各個點之間的最短路徑，A、（12）、B、C對應1到15這15個網(wǎng)點。例如要找A點到點的最短路，就是從path矩陣尋找。A到，即為1到8，首先找到矩陣中點（1，8）為數(shù)字12；再從12出發(fā)，找到點（12，8）數(shù)字為6；再從6出發(fā)，找到點（6，8）為數(shù)字15；最后從15出發(fā)，找到點（15，8）為數(shù)字8。所以最優(yōu)路線即為從1-12-

19、6-15-8，即從A到（11），（11）到，到C，再從C到，是從A點到點的最短路，其長度可以直接從圖3中按（1,8）找得，為22。于是很用以得到圖3中紅色框內(nèi)的部分分別就是從A、B、C到各菜市場的最短距離，整理出表3所示的每一百公斤蔬菜從各收購點到各菜市場的運輸費用。表3 每一百公斤蔬菜從各收購點到各菜市場的運輸費用（元）點購收場市菜A488191162220B14771612162317C20191114615510由于每天三個收購點的總供應量為200+170+160=530（100kg）每天8個菜市場的總需求量為75+60+80+70+100+55+90+80=610（100kg）所以每天

20、的短缺損失量為610-530=80（100kg）（一）對于（a）問題，可以建立以下lindo程序下的線性規(guī)劃模型：min 4a1+8b1+8c1+19d1+11e1+6f1+22g1+20h1+14a2+7b2+7c2+16d2+12e2+16f2+23g2+17h2+20a3+19b3+11c3+14d3+6e3+15f3+5g3+10h3+10a+8b+5c+10d+10e+8f+5g+8hst 1) a1+b1+c1+d1+e1+f1+g1+h1=2002)a2+b2+c2+d2+e2+f2+g2+h2=1703)a3+b3+c3+d3+e3+f3+g3+h3=1604)a+b+c+d+

21、e+f+g+h=805)a1+a2+a3+a=756)b1+b2+b3+b=607)c1+c2+c3+c=808)d1+d2+d3+d=709)e1+e2+e3+e=10010)f1+f2+f3+f=5511)g1+g2+g3+g=9012)h1+h2+h3+h=80End根據(jù)附錄1里面的運行結果，我們可以得出各收購點到各菜市場的蔬菜調運量，表4，其最小的蔬菜調運費用及預期短缺損失共為4610元。表4 蔬菜調運量（百公斤）點購收場市菜A750400305500B06040700000C0000700900短缺000000080由于有些收購點到菜市場的最短距離不是直接到達，而是經(jīng)過其他中轉站、菜

22、市場，甚至其他收購點的，因此還要根據(jù)圖4查出具體的調運線路如下：P（A，1）=A-1，運送量為75百公斤；P（A，3）=A-3，運送量為40百公斤；P（A，5）=A-(11)-5，運送量為30百公斤；P（A，6）=A-6，運送量為55百公斤；P（B，2）=B-2，運送量為60百公斤；P（B，3）=B-3，運送量為40百公斤；P（B，4）=B-（12）-4，運送量為70百公斤；P（C，5）=C-5，運送量為70百公斤；P（C，7）=C-7，運送量為90百公斤；注：P（i,j）為i到j的最短路徑因此，按點和點來表示，調運方案還可以這樣：表5 點對點的調運方案起點終點（11）（12）A7540553

23、0B604070C7090（11）30（12）70表中單位均為（百公斤），（11）、（12）為中轉點，空白處為零，表示無運輸需要。（二）顯然上一模型得出的結果中菜市場完全沒有得到供應，現(xiàn)實生活中是不允許的，那么接下來就解答第二個問題：若規(guī)定各菜市場短缺量一律不超過需求量的20%，重新設計定點供應方案。同樣的可以用線性模型來求解，于是建立以下lindo程序下的線性規(guī)劃模型：min 4a1+8b1+8c1+19d1+11e1+6f1+22g1+20h1+14a2+7b2+7c2+16d2+12e2+16f2+23g2+17h2+20a3+19b3+11c3+14d3+6e3+15f3+5g3+10

24、h3+10a+8b+5c+10d+10e+8f+5g+8hst 1) a1+b1+c1+d1+e1+f1+g1+h1=2002)a2+b2+c2+d2+e2+f2+g2+h2=1703)a3+b3+c3+d3+e3+f3+g3+h3=1604)a+b+c+d+e+f+g+h=805)a1+a2+a3+a=756)b1+b2+b3+b=607)c1+c2+c3+c=808)d1+d2+d3+d=709)e1+e2+e3+e=10010)f1+f2+f3+f=5511)g1+g2+g3+g=9012)h1+h2+h3+h=8013)a<1514)b<1215)c<1616)d&l

25、t;1417)e<2018)f<1119)g<1820)h<16End根據(jù)附錄2里面的運行結果，我們可以得出各收購點到各菜市場的蔬菜調運量，表6，其最小的蔬菜調運費用及預期短缺損失共為4806元。表6 蔬菜調運量（百公斤）點購收場市菜A751000606500B05064560000C00002407264短缺0016141601816同樣，根據(jù)圖4查出具體的調運線路如下：P（A，1）=A-1，運送量為75百公斤；P（A，2）=A-2，運送量為10百公斤；P（A，5）=A-(11)-5，運送量為60百公斤；P（A，6）=A-6，運送量為65百公斤；P（B，2）=B-2，

26、運送量為50百公斤；P（B，3）=B-3，運送量為64百公斤；P（B，4）=B-（12）-4，運送量為56百公斤；P（C，5）=C-5，運送量為24百公斤；P（C，7）=C-7，運送量為72百公斤；P（C，8）=C-8，運送量為64百公斤。因此，按點和點來表示，調運方案如下表7 點對點的調運方案（百公斤）起點終點1112A75100006500600B0506400000056C000024072640011000060000001200056000000（三）要滿足城市居民的蔬菜供應同時又要符合經(jīng)濟，就必須是理想的收購量總和等于需求總和，從表7得知，目前的最優(yōu)運輸線路中沒有一個收購點是同時作

27、為其他運輸線的中轉站的，因此，只需把新增加的蔬菜按最優(yōu)方案加到原來的基礎上，而不需要把原來的收購量做減少，可以建立以下lindo程序下的線性規(guī)劃模型：min 4a1+8b1+8c1+19d1+11e1+6f1+22g1+20h1+14a2+7b2+7c2+16d2+12e2+16f2+23g2+17h2+20a3+19b3+11c3+14d3+6e3+15f3+5g3+10h3st1)a1+a2+a3=752)b1+b2+b3=603)c1+c2+c3=804)d1+d2+d3=705)e1+e2+e3=1006)f1+f2+f3=557)g1+g2+g3=908)h1+h2+h3=809)

28、a1+b1+c1+d1+e1+f1+g1+h1>20010)a2+b2+c2+d2+e2+f2+g2+h2>17011)a3+b3+c3+d3+e3+f3+g3+h3>160end從附錄3的運算結果可以得知，增加的80百公斤手工量只需全部分給C收購點，然后重新分配，如表8所示，這時滿足了題目要求，同時，蔬菜調運費用及預期短缺損失也是最低的，共4770元。表8 蔬菜調運量（百公斤）點購收場市菜A754000305500B02080700000C00007009080短缺00000000六、模型的檢驗和改進本文所涉及的最短路問題以及運輸供求平衡問題都沒有使用常規(guī)的圖上作業(yè)法，而是

29、把算法的思路轉化為計算機程序，節(jié)省了求解的時間同時也提高了模型的準確度，而且具有較強的可復制性。當然在設計最少運費（最短路）的時候，可以嘗試把它和接下來的運輸問題用“供求平衡”的方法結合在一起求解，即不用先求出各采購點到菜市場的最短路徑，直接按必須供應部分和選擇性供應部分用動態(tài)規(guī)劃的方法求解以驗證，由于手工的運算量大，在此不作贅述。而第三問的增加供應分配方案，可以從附錄2的靈敏度分析中得知改變?nèi)魏我粋€的供應量都接改變最終的結果（因為ROW1、2、3的ALLOWABLE INCREASE和ALLOWABLE DECREASE都為零），研究影子價格的意義就不大了。綜上所述，本文中所涉及的方法和模型

30、都是合適的。七、模型的推廣 本文的解題思路以及所涉及的方法和模型都是準確而且可復制性強的，在解決各種最小費用、最短路線、產(chǎn)銷平衡、運輸問題時都有較強的參考意義，適當?shù)倪\用計算機程序解決復雜的計算問題有利于數(shù)學使用的推廣。八、參考文獻【1】運籌學教材編寫組，運籌學，清華大學出版社，2011九、附錄1、lindo運算結果1（含靈敏度分析） LP OPTIMUM FOUND AT STEP 1 OBJECTIVE FUNCTION VALUE 1) 4610.000 VARIABLE VALUE REDUCED COST A1 75.000000 0.000000 B1 0.000000 0.000

31、000 C1 40.000000 0.000000 D1 0.000000 2.000000 E1 30.000000 0.000000 F1 55.000000 0.000000 G1 0.000000 12.000000 H1 0.000000 5.000000 A2 0.000000 11.000000 B2 60.000000 0.000000 C2 40.000000 0.000000 D2 70.000000 0.000000 E2 0.000000 2.000000 F2 0.000000 11.000000 G2 0.000000 14.000000 H2 0.000000 3

32、.000000 A3 0.000000 21.000000 B3 0.000000 16.000000 C3 0.000000 8.000000 D3 0.000000 2.000000 E3 70.000000 0.000000 F3 0.000000 14.000000 G3 90.000000 0.000000 H3 0.000000 0.000000 A 0.000000 13.000000 B 0.000000 7.000000 C 0.000000 4.000000 D 0.000000 0.000000 E 0.000000 6.000000 F 0.000000 9.00000

33、0 G 0.000000 2.000000 H 80.000000 0.000000 ROW SLACK OR SURPLUS DUAL PRICES 1) 0.000000 -15.000000 2) 0.000000 -14.000000 3) 0.000000 -10.000000 4) 0.000000 -8.000000 5) 0.000000 11.000000 6) 0.000000 7.000000 7) 0.000000 7.000000 8) 0.000000 -2.000000 9) 0.000000 4.000000 10) 0.000000 9.000000 11)

34、0.000000 5.000000 12) 0.000000 0.000000 NO. ITERATIONS= 1 RANGES IN WHICH THE BASIS IS UNCHANGED: OBJ COEFFICIENT RANGES VARIABLE CURRENT ALLOWABLE ALLOWABLE COEF INCREASE DECREASE A1 4.000000 11.000000 INFINITY B1 8.000000 INFINITY 0.000000 C1 8.000000 0.000000 2.000000 D1 19.000000 INFINITY 2.0000

35、00 E1 11.000000 2.000000 0.000000 F1 6.000000 9.000000 INFINITY G1 22.000000 INFINITY 12.000000 H1 20.000000 INFINITY 5.000000 A2 14.000000 INFINITY 11.000000 B2 7.000000 0.000000 INFINITY C2 7.000000 2.000000 0.000000 D2 16.000000 0.000000 2.000000 E2 12.000000 INFINITY 2.000000 F2 16.000000 INFINI

36、TY 11.000000 G2 23.000000 INFINITY 14.000000 H2 17.000000 INFINITY 3.000000 A3 20.000000 INFINITY 21.000000 B3 19.000000 INFINITY 16.000000 C3 11.000000 INFINITY 8.000000 D3 14.000000 INFINITY 2.000000 E3 6.000000 0.000000 2.000000 F3 15.000000 INFINITY 14.000000 G3 5.000000 2.000000 INFINITY H3 10.

37、000000 INFINITY 0.000000 A 10.000000 INFINITY 13.000000 B 8.000000 INFINITY 7.000000 C 5.000000 INFINITY 4.000000 D 10.000000 2.000000 0.000000 E 10.000000 INFINITY 6.000000 F 8.000000 INFINITY 9.000000 G 5.000000 INFINITY 2.000000 H 8.000000 0.000000 INFINITY RIGHTHAND SIDE RANGES ROW CURRENT ALLOW

38、ABLE ALLOWABLE RHS INCREASE DECREASE 1 200.000000 0.000000 0.000000 2 170.000000 0.000000 0.000000 3 160.000000 0.000000 0.000000 4 80.000000 0.000000 0.000000 5 75.000000 0.000000 0.000000 6 60.000000 0.000000 0.000000 7 80.000000 0.000000 0.000000 8 70.000000 0.000000 0.000000 9 100.000000 0.00000

39、0 0.000000 10 55.000000 0.000000 0.000000 11 90.000000 0.000000 0.000000 12 80.000000 0.000000 0.0000002、lindo運算結果2（含靈敏度分析） LP OPTIMUM FOUND AT STEP 20 OBJECTIVE FUNCTION VALUE 1) 4806.000 VARIABLE VALUE REDUCED COST A1 75.000000 0.000000 B1 10.000000 0.000000 C1 0.000000 0.000000 D1 0.000000 2.000000 E1 60.000000 0.000000 F1 55.000000 0.000000 G1 0.000000 12.000000 H1 0.000000 5.000000 A2 0.000000 11.000000 B2 50.000000 0.000000 C2 64.000000 0.000000 D2 56.000000 0.000000 E2 0.000000 2.000000 F2 0.000000 11.000000 G2 0.000000 1