?前面介紹了利用Matlab二元整數規劃求解數獨問題，對于另一個問題-旅行商問題也可以用它來求解。

旅行商問題就是找到經過所有站點的最短閉合路徑，如下圖為在美國地圖框架內產生的200個旅行站點，而旅行商要找到一條最短路徑將200個站點都旅行到。這也可以借助二元整數規劃求解。 這個例子比較典型： 首先，例子的計算規模很大，200個站點，相應就有19900個二元變量，對如此多的變量同時進行優化，計算量可想而知，所以求解中使用了稀疏矩陣表達約束，也就是說Matlab優化函數是接受稀疏矩陣形式的約束條件的； ? ??另外，對于常規的求解過程，會出現局部旅行問題，如下圖出現的獨立閉合圈。為了避免這些閉合圈的出現，通常通過增加不等式約束的方式解決，但相關不等式約束的數量將非常龐大，因為任意幾個點之間都要保證不出現閉合。為此，例子使用了不斷增加不等式約束的迭代優化方式，即先不添加不等式約束進行優化，如果優化結果中出現閉合圈，則添加相應的約束條件避免這些圈的存在，然后再進行優化，直到不再有閉合圈存在。
? ? 這個例子包含3個文件：TravellingSalesmanExample.m、updateSalesmanPlot.m和detectSubtours.m。其中第一個是整個問題的求解過程，可以用于發布，第二個用于更新每步的優化結果圖，第三個為識別每次結果的閉合圈。下面將注釋翻譯后的文件粘貼如下 ? ??TravellingSalesmanExample.m文件： %%?旅行商問題，Travelling?Salesman?Problem %?這個例子展示如何使用二元整數規劃求解經典的旅行商問題。這個問題就是找到經過所有 %?站點的最短閉合路徑。這個例子使用了200個站點，但也可以通過改變nStops來得到不同 %?運算量的求解過程。 %?對于初始情況，求解結果中存在子旅途，即結果沒有給出一條經過所有站點的連續路線， %?而是有一些獨立的環存在。對此，可以通過迭代加約束、再優化的方法實現子旅途的去除。 ? clc clear close?all ? %%?繪制地圖和站點 nStops?=?200;?%?站點個數，可以自己調整，不過求解復雜度與其平方成正比 stopsLon?=?zeros(nStops,1);?%?分配站點的x坐標 stopsLat?=?stopsLon;?%?分配y坐標 ? %?在代表美國的粗糙多邊形內產生隨機站點 load('usborder.mat','x','y','xx','yy');?%?加載數據，其中(x,y)為美國的精細輪廓，(xx,yy)為粗糙輪廓 rng(3,'twister')?%?初始化隨機序列，保證每次運行的重復性 n?=?1; while?(n?<=?nStops) ????xp?=?rand*1.5; ????yp?=?rand; ???? ????%?是否在(xx,yy)的粗糙多邊形內 ????if?inpolygon(xp,yp,xx,yy) ????????stopsLon(n)?=?xp; ????????stopsLat(n)?=?yp; ????????n?=?n+1; ????end end ? figure plot(x,y,'Color','red');?%?多邊形邊界 hold?on plot(stopsLon,stopsLat,'*b')?%?站點 hold?off pause(1) ? %%?問題建模 %?將旅行商問題建模為整數線性規劃問題： %?1.?得到所有可能的旅行，即所有站點對； %?2.?計算各旅行的距離； %?3.?代價函數就是最小化旅途距離； %?4.?決定變量是二元的，即如果旅行出現在旅途中，則為1，否則為0； %?5.?為了保證旅途包含所有站點，可以設定如下線性約束：每個站點出現在兩個旅行中，即達到和離開。 ? %%?計算點之間的距離 %?如果有200個站點，就有19900個旅行，即19900個二元變量。 idxs?=?nchoosek(1:nStops,2);?%?得到所有組合 ? %%?計算旅行的距離 %?假設地球是平的，這樣可以快速計算 dist?=?hypot(stopsLat(idxs(:,1))?-?stopsLat(idxs(:,2)),?... ????stopsLon(idxs(:,1))?-?stopsLon(idxs(:,2))); lendist?=?length(dist); ? %%?小結1 %?基于上面定義的dist向量，旅途長度為dist'*x，其中x為二元解向量。這正是我們要最小化的代價函數。 ? %%?等式約束 %?問題有兩類等式約束，第一個保證共有200個旅行，第二個保證每個站點必須有2個旅行與其相關。 %?第一個可以用Aeq*x?=?beq的形式指定。 Aeq?=?spones(1:length(idxs));?%?全1向量，即將所有x相加 beq?=?nStops; ? %?第二個通過擴展上面的Aeq來完成 Aeq?=?[Aeq;spalloc(nStops,length(idxs),nStops*(nStops-1))];?%?分配一個稀疏矩陣 for?ii?=?1:nStops ????whichIdxs?=?(idxs?==?ii);?%?找到包含站點ii的所有旅行 ????whichIdxs?=?sparse(sum(whichIdxs,2));?%?只要有任何一端為ii，即包括此旅行 ????Aeq(ii+1,:)?=?whichIdxs';?%?加入Aeq矩陣 end beq?=?[beq;?2*ones(nStops,1)]; ? %%?二元邊界 %?所有變量都是二元的，所以設置intcon保證所有變量為整數，同時設置下邊界為0，上邊界為1 intcon?=?1:lendist; lb?=?zeros(lendist,1); ub?=?ones(lendist,1); ? %%?使用intlinprog完成最優化 opts?=?optimoptions('intlinprog','Display','off'); tic [xopt,costopt,exitflag,output]?=?intlinprog(dist,intcon,[],[],Aeq,beq,lb,ub,opts); toc ? %%?繪制結果 hold?on lh?=?zeros(nStops,1);?%?用來保存線的句柄 lh?=?updateSalesmanPlot(lh,xopt,idxs,stopsLon,stopsLat); title('帶有子旅途的結果'); ? %%?小結2 %?從繪制結果可以看到，圖中存在許多子旅途。到目前還沒有設置約束阻止這些子旅途的出現。 %?而為了阻止它們的出現，需要設置難以想象的大量不等式約束。 ? %%?子旅途約束 %?由于不可能一下添加所有子旅途約束，所以可以通過迭代的方式添加約束。 %?即在當前結果中發現子旅途時，則添加相應的不等式約束來阻止這些特定子旅途的出現。 %?通過這種方式，通過簡單的幾次迭代就可找到合適的旅途。 ? %?利用不等式約束消除子旅途，比如出現一個5點的子旅途，則有5條線將這些點連起來， %?所以可以要求這五個點之間的線不多于4條。如果找到這五個點之間的所有線，然后限定解中最多出現4條這種線。 %?這樣限定是合理的，因為如果有5個或更多線出現在結果中，則將必然出現子旅途(具有n個節點和邊的圖必然存在圈)。 ? %?函數detectSubtours通過分析求解結果返回向量元胞數組，每個向量包含了組成子旅途的站點。 tours?=?detectSubtours(xopt,idxs); numtours?=?length(tours);?%?子旅途的個數 fprintf('子旅途的數目為:?%d\n',?numtours); ? %%?加入線性不等式約束消除子旅途，然后重復優化過程，直到只得到一個子旅途為止 A?=?spalloc(0,lendist,0);?%?分配稀疏線性不等式約束矩陣 b?=?[]; rowIdx?=?0; while?numtours?>?1?%?循環直到只有一個子旅途 ????%?添加子旅途約束 ????b?=?[b;zeros(numtours,1)]; ????A?=?[A;spalloc(numtours,lendist,nStops)];?%?初始分配的稀疏空間不那么重要，這里設為nStops，少了會自動分配空間 ????for?ii?=?1:numtours?%?添加numtours個不等式約束 ????????rowIdx?=?rowIdx+1;?%?指向約束的指針 ????????subTourIdx?=?tours{ii};?%?提取當前子旅途 ? ????????%?找到與當前子旅途關聯的所有變量，然后添加不等式約束來阻止所有途徑這些站點的子旅途 ????????variations?=?nchoosek(1:length(subTourIdx),2); ????????for?jj?=?1:length(variations) ????????????whichVar?=?(sum(idxs==subTourIdx(variations(jj,1)),2))?&?... ????????????????(sum(idxs==subTourIdx(variations(jj,2)),2)); ????????????A(rowIdx,whichVar)?=?1; ????????end ????????b(rowIdx)?=?length(subTourIdx)-1;?%?子旅途站點減1為b的值 ????end ???? ????%?再次嘗試優化 ????tic ????[xopt,costopt,exitflag,output]?=?intlinprog(dist,intcon,A,b,Aeq,beq,lb,ub,opts); ????toc ???? ????%?顯示結果 ????lh?=?updateSalesmanPlot(lh,xopt,idxs,stopsLon,stopsLat); ???? ????%?檢測這次的子旅途數目? ????tours?=?detectSubtours(xopt,idxs); ????numtours?=?length(tours); ????fprintf('子旅途的數目:?%d\n',numtours); end ? title('除去子旅途的結果'); hold?off ? %%?求解效果 %?結果給出了一個可行旅途，但它是否是最小代價的旅途呢？一種檢驗方式就是檢查輸出變量output output? %?output中顯示內部計算的上下限間隔為0，所以結果為最優解。? ? ? ??updateSalesmanPlot.m文件： function?lh?=?updateSalesmanPlot(lh,xopt,idxs,stopsLat,stopsLon) %?tsp_intlinprog例子的繪圖函數 ? if?(?lh?~=?zeros(size(lh))?)?%?通過lh是否全零來判斷是否第一次進入 ????set(lh,'Visible','off');?%?移除之前的線 end ? %?求解結果中的旅行 segments?=?find(xopt);? ? %?形成繪制數據 Lat?=?zeros(3*length(segments),1); Lon?=?zeros(3*length(segments),1); for?ii?=?1:length(segments) ????start?=?idxs(segments(ii),1); ????stop?=?idxs(segments(ii),2); ???? ????%?數據之間用NaN分割，這樣可以繪制分離的線段，不然需要一個一個繪制線段，產生很多句柄 ????Lat(3*ii-2:3*ii)?=?[stopsLat(start);?stopsLat(stop);?NaN]; ????Lon(3*ii-2:3*ii)?=?[stopsLon(start);?stopsLon(stop);?NaN]; end ? lh?=?plot(Lat,Lon,'k:','LineWidth',2); set(lh,'Visible','on'); drawnow; ? ? ??detectSubtours.m文件： function?subTours?=?detectSubtours(x,idxs) %?返回子旅途元胞數組，即圖中的子圈 ? x?=?round(x);?%?糾正不確切的整數 r?=?find(x); substuff?=?idxs(r,:);?%?旅行線的節點對 unvisited?=?ones(length(r),1);?%?跟蹤未訪問的旅行 curr?=?1;?%?正在評價的子旅途 startour?=?find(unvisited,1);?%?第一個未訪問的旅行 while?~isempty(startour) ????home?=?substuff(startour,1); ????nextpt?=?substuff(startour,2); ????visited?=?nextpt; ????unvisited(startour)?=?0; ???? ????while?nextpt?~=?home ????????%?找以nextpt為起點的旅行 ????????[srow,scol]?=?find(substuff?==?nextpt); ???????? ????????%?確定相應旅行的節點 ????????trow?=?srow(srow?~=?startour); ????????scol?=?3-scol(trow?==?srow);?%?1變2，2變1 ????????startour?=?trow; ????????nextpt?=?substuff(startour,scol);?%?子旅途的下一節點位置 ???????? ????????visited?=?[visited,nextpt];?%?將節點加入子旅途 ????????unvisited(startour)?=?0;?%?更新訪問過的位置 ????end ????subTours{curr}?=?visited;?%?保存找到的子旅途 ???? ????curr?=?curr?+?1; ????startour?=?find(unvisited,1); end end ? 運行TravellingSalesmanExample.m，得到的結果如下 Elapsed time is 3.394080 seconds. 子旅途的數目為: 28 Elapsed time is 12.648189 seconds. 子旅途的數目: 22 Elapsed time is 1.452217 seconds. 子旅途的數目: 10 Elapsed time is 3.314973 seconds. 子旅途的數目: 5 Elapsed time is 4.305792 seconds. 子旅途的數目: 3 Elapsed time is 8.356620 seconds. 子旅途的數目: 4 Elapsed time is 16.728625 seconds. 子旅途的數目: 1 output = relativegap: 0.0100 absolutegap: 0.0021 numfeaspoints: 5 numnodes: 2258 constrviolation: 2.5580e-13 message: 'Optimal solution found. Intlinprog stopped because the objective value is wi...'
from:?http://chunqiu.blog.ustc.edu.cn/?p=812

總結

以上是生活随笔為你收集整理的利用Matlab优化工具箱求解旅行商最短路径问题的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。