基于改进遗传算法求解 TSP/MTSP 【旅行推销员问题（TSP）、多旅行推销员问题（M-TSP）】（Matlab代码实现）

1 概述

2 旅行商问题

3 遗传算法

4 约束优化

5 带有罚方法的遗传算法的流程图

6 带有惩罚函数的遗传算法在TSP中的应用

7 运行结果

7.1 单旅行商问题

7.2 多旅行商问题

8 Matlab代码实现

9 参考文献

10 写在最后

1 概述

主要研究用遗传算法解决带有约束的TSP的方法。使用贪婪交叉算子、自适应变异算子和带有精英保留策略的选择算子相结合对基本遗传算法进行了改进,针对实际TSP中的约束条件讨论了罚方法在遗传算法中的应用,提出了自适应的惩罚函数,并将其与改进后的遗传算法相结合,解决了带有时间约束的TSP。通过对实验结果的比较分析,证明了该方法的可行性和有效性。

2 旅行商问题

3 遗传算法

遗传算法是1975年由John Holland教授提出。它借用了仿真生物遗传学和自然选择机理,通过自然选择、交叉、变异等机制实现各个个体适应性的提高。遗传算法在整个进化过程中的遗传操作是随机的;但是它呈现出来的特性却并不是完全随机的,它能够有效地利用历史信息使得下一代群体性能有所提高,一代一代不断进化,最终收敛到一个即使不是最优解,也是离最优解非常接近的解空间上。遗传算法涉及到编码方式、初始种群、适应度函数、遗传算子和控制参数五大要素。
基本遗传算法虽然具有全局搜索能力,但是在进化过程中往往容易出现收敛速度缓慢或不收敛,又或陷入局部最优的情况,加快收敛速度又有可能出现早熟现象,过或不及都不能及时准确地获得最优解。因此,学者们在这些方面进行了大量的研究,希望既能加快收敛速度,又能防止早熟收敛。

4 约束优化

5 带有罚方法的遗传算法的流程图

6 带有惩罚函数的遗传算法在TSP中的应用

取一段路线并颠倒该段中城市的顺序

取路线中的两个城市并按路径顺序交换它们的位置

将一个城市拉出路线并将其重新插入另一个地方，将其他城市滑落

7 运行结果

7.1 单旅行商问题

7.2 多旅行商问题

8 Matlab代码实现

本文仅展现部分代码，全部代码点击链接：

基于改进遗传算法求解 TSP/MTSP 【旅行推销员问题（TSP）、多旅行推销员问题（M-TSP）】（Matlab代码实现）

function varargout = mtspof_ga(varargin)%% Initialize default configuration%defaultConfig.xy          = 10*rand(40,2);defaultConfig.dmat        = [];defaultConfig.nSalesmen   = 5;defaultConfig.minTour     = 1;defaultConfig.popSize     = 80;defaultConfig.numIter     = 5e3;defaultConfig.showProg    = true;defaultConfig.showStatus  = true;defaultConfig.showResult  = true;defaultConfig.showWaitbar = false;%% Interpret user configuration inputs%if ~narginuserConfig = struct();elseif isstruct(varargin{1})userConfig = varargin{1};elsetryuserConfig = struct(varargin{:});catcherror('??? Expected inputs are either a structure or parameter/value pairs');endend%% Override default configuration with user inputs%configStruct = get_config(defaultConfig,userConfig);%% Extract configuration%xy          = configStruct.xy;dmat        = configStruct.dmat;nSalesmen   = configStruct.nSalesmen;minTour     = configStruct.minTour;popSize     = configStruct.popSize;numIter     = configStruct.numIter;showProg    = configStruct.showProg;showStatus  = configStruct.showStatus;showResult  = configStruct.showResult;showWaitbar = configStruct.showWaitbar;if isempty(dmat)nPoints = size(xy,1);a = meshgrid(1:nPoints);dmat = reshape(sqrt(sum((xy(a,:)-xy(a',:)).^2,2)),nPoints,nPoints);end%% Verify inputs%[N,dims] = size(xy);[nr,nc] = size(dmat);if (N ~= nr) || (N ~= nc)error('??? Invalid XY or DMAT inputs')endn = N - 2; % Separate start and end cities%% Sanity checks%nSalesmen   = max(1,min(n,round(real(nSalesmen(1)))));minTour     = max(1,min(floor(n/nSalesmen),round(real(minTour(1)))));popSize     = max(8,8*ceil(popSize(1)/8));numIter     = max(1,round(real(numIter(1))));showProg    = logical(showProg(1));showStatus  = logical(showStatus(1));showResult  = logical(showResult(1));showWaitbar = logical(showWaitbar(1));%% Initializations for route break point selection%nBreaks = nSalesmen-1;dof = n - minTour*nSalesmen;          % degrees of freedomaddto = ones(1,dof+1);for k = 2:nBreaksaddto = cumsum(addto);endcumProb = cumsum(addto)/sum(addto);%% Initialize the populations%popRoute = zeros(popSize,n);         % population of routespopBreak = zeros(popSize,nBreaks);   % population of breakspopRoute(1,:) = (1:n) + 1;popBreak(1,:) = rand_breaks();for k = 2:popSizepopRoute(k,:) = randperm(n) + 1;popBreak(k,:) = rand_breaks();end%% Seed the algorithm with a previous result if available%if all(isfield(userConfig,{'optRoute','optBreak'}))optRoute = userConfig.optRoute;optBreak = userConfig.optBreak;isValidRoute = isequal(popRoute(1,:),sort(optRoute));isValidBreak = all(optBreak > 0) && all(optBreak <= n) && ...(length(optBreak) == nBreaks) && ~any(mod(optBreak,1));if isValidRoute && isValidBreakpopRoute(1,:) = optRoute;popBreak(1,:) = optBreak;endend%% Select the colors for the plotted routes%pclr = ~get(0,'DefaultAxesColor');clr = [1 0 0; 0 0 1; 0.67 0 1; 0 1 0; 1 0.5 0];if (nSalesmen > 5)clr = hsv(nSalesmen);end%% Run the GA%row = zeros(popSize,n+nSalesmen);col = zeros(popSize,n+nSalesmen);isValid = false(1,n+nSalesmen);globalMin = Inf;distHistory = NaN(1,numIter);tmpPopRoute = zeros(8,n);tmpPopBreak = zeros(8,nBreaks);newPopRoute = zeros(popSize,n);newPopBreak = zeros(popSize,nBreaks);[isClosed,isStopped,isCancelled] = deal(false);if showProghFig = figure('Name','MTSPOF_GA | Current Best Solution', ...'Numbertitle','off','CloseRequestFcn',@close_request);hAx = gca;if showStatus[hStatus,isCancelled] = figstatus(0,numIter,[],hFig);endendif showWaitbarhWait = waitbar(0,'Searching for near-optimal solution ...', ...'CreateCancelBtn',@cancel_search);endisRunning = true;for iter = 1:numIter%% EVALUATE SOLUTIONS%   This section of code computes the total cost of each solution%   in the population. The actual code that gets executed uses a%   much faster (vectorized) method to calculate the route lengths%   compared to the triple for-loop below (provided for reference)%   but gives the same result.%%     totalDist = zeros(popSize,1);%     for p = 1:popSize%         d = 0;%         pRoute = popRoute(p,:);%         pBreak = popBreak(p,:);%         rng = [[1 pBreak+1];[pBreak n]]';%         for s = 1:nSalesmen%             d = d + dmat(1,pRoute(rng(s,1)));%             for k = rng(s,1):rng(s,2)-1%                 d = d + dmat(pRoute(k),pRoute(k+1));%             end%             d = d + dmat(pRoute(rng(s,2)),N);%         end%         totalDist(p) = d;%     end%for p = 1:popSizebrk = popBreak(p,:);isValid(:) = false;isValid([1 brk+(2:nSalesmen)]) = true;row(p,isValid) = 1;row(p,~isValid) = popRoute(p,:);isValid(:) = false;isValid([brk+(1:nSalesmen-1) n+nSalesmen]) = true;col(p,isValid) = N;col(p,~isValid) = popRoute(p,:);endind = N*(col-1) + row;totalDist = sum(dmat(ind),2);%% SELECT THE BEST%   This section of code finds the best solution in the current%   population and stores it if it is better than the previous best.%[minDist,index] = min(totalDist);distHistory(iter) = minDist;if (minDist < globalMin)globalMin = minDist;optRoute = popRoute(index,:);optBreak = popBreak(index,:);rng = [[1 optBreak+1];[optBreak n]]';if showProg%% Plot the best route%for s = 1:nSalesmenrte = [1 optRoute(rng(s,1):rng(s,2)) N];if (dims > 2), plot3(hAx,xy(rte,1),xy(rte,2),xy(rte,3),'.-','Color',clr(s,:));else, plot(hAx,xy(rte,1),xy(rte,2),'.-','Color',clr(s,:)); endhold(hAx,'on');endif (dims > 2), plot3(hAx,xy(1,1),xy(1,2),xy(1,3),'o','Color',pclr);plot3(hAx,xy(N,1),xy(N,2),xy(N,3),'o','Color',pclr);else, plot(hAx,xy(1,1),xy(1,2),'o','Color',pclr);plot(hAx,xy(N,1),xy(N,2),'o','Color',pclr);endtitle(hAx,sprintf('Total Distance = %1.4f, Iteration = %d',minDist,iter));hold(hAx,'off');drawnow;endend%% Update the status bar and check cancellation status%if showProg && showStatus && ~mod(iter,ceil(numIter/100))[hStatus,isCancelled] = figstatus(iter,numIter,hStatus,hFig);endif (isStopped || isCancelled)breakend%% MODIFY THE POPULATION%   This section of code invokes the genetic algorithm operators.%   In this implementation, solutions are randomly assigned to groups%   of eight and the best solution is kept (tournament selection).%   The best-of-eight solution is then mutated 3 different ways%   (flip, swap, and slide) and the four resulting solutions are%   then copied but assigned different break points to complete the%   group of eight. There is no crossover operator because it tends%   to be highly destructive and rarely improves a decent solution.%randomOrder = randperm(popSize);for p = 8:8:popSizertes = popRoute(randomOrder(p-7:p),:);brks = popBreak(randomOrder(p-7:p),:);dists = totalDist(randomOrder(p-7:p));[ignore,idx] = min(dists); %#okbestOf8Route = rtes(idx,:);bestOf8Break = brks(idx,:);routeInsertionPoints = sort(randperm(n,2));I = routeInsertionPoints(1);J = routeInsertionPoints(2);for k = 1:8 % Generate new solutionstmpPopRoute(k,:) = bestOf8Route;tmpPopBreak(k,:) = bestOf8Break;switch kcase 2 % FliptmpPopRoute(k,I:J) = tmpPopRoute(k,J:-1:I);case 3 % SwaptmpPopRoute(k,[I J]) = tmpPopRoute(k,[J I]);case 4 % SlidetmpPopRoute(k,I:J) = tmpPopRoute(k,[I+1:J I]);case 5 % Modify breakstmpPopBreak(k,:) = rand_breaks();case 6 % Flip, modify breakstmpPopRoute(k,I:J) = tmpPopRoute(k,J:-1:I);tmpPopBreak(k,:) = rand_breaks();case 7 % Swap, modify breakstmpPopRoute(k,[I J]) = tmpPopRoute(k,[J I]);tmpPopBreak(k,:) = rand_breaks();case 8 % Slide, modify breakstmpPopRoute(k,I:J) = tmpPopRoute(k,[I+1:J I]);tmpPopBreak(k,:) = rand_breaks();otherwise % Do nothingendendnewPopRoute(p-7:p,:) = tmpPopRoute;newPopBreak(p-7:p,:) = tmpPopBreak;endpopRoute = newPopRoute;popBreak = newPopBreak;%% Update the waitbar%if showWaitbar && ~mod(iter,ceil(numIter/325))waitbar(iter/numIter,hWait);endendif showProg && showStatusfigstatus(numIter,numIter,hStatus,hFig);endif showWaitbardelete(hWait);endisRunning = false;if isCloseddelete(hFig);end%% Append prior distance history if present%if isfield(userConfig,'distHistory')priorHistory = userConfig.distHistory;isNan = isnan(priorHistory);distHistory = [priorHistory(~isNan) distHistory];end%% Format the optimal solution%optSolution = cell(nSalesmen,1);rng = [[1 optBreak+1];[optBreak n]]';for s = 1:nSalesmenoptSolution{s} = [1 optRoute(rng(s,1):rng(s,2)) N];end%% Show the final results%if showResult%% Plot the GA results%figure('Name','MTSPOF_GA | Results','Numbertitle','off');subplot(2,2,1);if (dims > 2), plot3(xy(:,1),xy(:,2),xy(:,3),'.','Color',pclr);else, plot(xy(:,1),xy(:,2),'.','Color',pclr); endtitle('City Locations');subplot(2,2,2);imagesc(dmat([1 optRoute N],[1 optRoute N]));title('Distance Matrix');subplot(2,2,3);for s = 1:nSalesmenrte = optSolution{s};if (dims > 2), plot3(xy(rte,1),xy(rte,2),xy(rte,3),'.-','Color',clr(s,:));else, plot(xy(rte,1),xy(rte,2),'.-','Color',clr(s,:)); endtitle(sprintf('Total Distance = %1.4f',minDist));hold on;endif (dims > 2), plot3(xy(1,1),xy(1,2),xy(1,3),'o','Color',pclr);plot3(xy(N,1),xy(N,2),xy(N,3),'o','Color',pclr);else, plot(xy(1,1),xy(1,2),'o','Color',pclr);plot(xy(N,1),xy(N,2),'o','Color',pclr);endsubplot(2,2,4);plot(distHistory,'b','LineWidth',2);title('Best Solution History');set(gca,'XLim',[1 length(distHistory)],'YLim',[0 1.1*max([1 distHistory])]);end%% Return output%if nargout%% Create anonymous functions for plot generation%plotPoints  = @(s)plot(s.xy(:,1),s.xy(:,2),'.','Color',~get(gca,'Color'));plotResult  = @(s)cellfun(@(s,i)plot(s.xy(i,1),s.xy(i,2),'.-', ...'Color',rand(1,3)),repmat({s},size(s.optSolution)),s.optSolution);plotHistory = @(s)plot(s.distHistory,'b-','LineWidth',2);plotMatrix  = @(s)imagesc(s.dmat(cat(2,s.optSolution{:}),cat(2,s.optSolution{:})));%% Save results in output structure%resultStruct = struct( ...'xy',          xy, ...'dmat',        dmat, ...'nSalesmen',   nSalesmen, ...'minTour',     minTour, ...'popSize',     popSize, ...'numIter',     numIter, ...'showProg',    showProg, ...'showResult',  showResult, ...'showWaitbar', showWaitbar, ...'optRoute',    optRoute, ...'optBreak',    optBreak, ...'optSolution', {optSolution}, ...'plotPoints',  plotPoints, ...'plotResult',  plotResult, ...'plotHistory', plotHistory, ...'plotMatrix',  plotMatrix, ...'distHistory', distHistory, ...'minDist',     minDist);varargout = {resultStruct};end%% Generate random set of break points%function breaks = rand_breaks()if (minTour == 1) % No constraints on breaksbreaks = sort(randperm(n-1,nBreaks));else % Force breaks to be at least the minimum tour lengthnAdjust = find(rand < cumProb,1)-1;spaces = randi(nBreaks,1,nAdjust);adjust = zeros(1,nBreaks);for kk = 1:nBreaksadjust(kk) = sum(spaces == kk);endbreaks = minTour*(1:nBreaks) + cumsum(adjust);endend%% Nested function to cancel search%function cancel_search(varargin)isStopped = true;end%% Nested function to close the figure window%function close_request(varargin)if isRunning[isClosed,isStopped] = deal(true);isRunning = false;elsedelete(hFig);endendend