ts_tsp.m

解决TSP问题的TS算法MATLAB实现

function [result, shunxu]= TS_TSP(N, inStart, IsBianYi)
% 用禁忌搜索算法解决旅行商问题(可以加入变异算子，效果更好)

figure(1), close;
figure(2), close;
figure(3), close;
figure(4), close;

displayflag= 1; % 控制是否显示结果信息
% IsBianYi= 0; % 是否选择变异算子

% 城市数目
% N= 29;

% 获得城市坐标
[X, Y, Nonfit, bestxu]= TSP_DATA(N);

if Nonfit==1
    rand('state', 0);
    X= rand(1, N);
    rand('state', 1);
    Y= rand(1, N);
end;

% 计算出城市之间的距离矩阵
for i= 1:N
    for j= 1:N
        dis(i, j)= sqrt((X(i)-X(j))^2+(Y(i)-Y(j))^2);
    end;
end;
% dissort= sort(dis);
% disMin= dissort(2, :)

if ~Nonfit
    if displayflag
        figure(1); plot(X(bestxu), Y(bestxu), 'b*-'); title('best order of cities')
    end;
    bestdis= 0;
    for i= 1:N
        bestdis= bestdis + dis(bestxu(i), bestxu(i+1));
    end;
    text(min(X)+0.9*(max(X)-min(X)), min(Y)+0.1*(max(Y)-min(Y)), num2str(bestdis));
end;

T= N*20;
TT= ceil(T/2);
tabulength= 5+N; % 禁忌长度 % 
runcount= 0; % 搜索次数 计数变量
Pby= 0.4;
Max= 100000;

%====
xulie= 1:N; % 总共的城市 初始序列
Start= inStart; % 起始城市
xulie(find(xulie==Start))= []; % 剩余城市
% % 随机排列剩余的城市的先后顺序
% % rand('state', 2); % 调试时使用
Maxeye= Max*eye(N);
distemp= dis + Maxeye;
distemp(:, Start)= Max;
xulie_1= find(distemp(Start, :)==min(distemp(Start, :)), 1, 'first');
distemp(:, xulie_1)= Max;
xulie_end= find(distemp(Start, :)==min(distemp(Start, :)), 1, 'first');
getrandom= rand(1, N-1-2);
sortrandom= sort(getrandom);
for randcnt= 1:N-1-2
    sortnum(randcnt)= find(getrandom==sortrandom(randcnt));
end;
xulie_temp= xulie;
xulie_temp([1, end])= [];
candidate_now.xulie= [xulie_1, xulie_temp(sortnum), xulie_end]; % 随机排序后的剩余城市  当前解（城市顺序）

%% 按照最近距离法则生成城市初始序列
% Maxeye= Max*eye(N);
% distemp= dis + Maxeye;
% distemp(:, Start)= Max;
% candidate_now.xulie(1)= find(distemp(Start, :)==min(distemp(Start, :)), 1, 'first');
% for i= 2:N-1
%     distemp(:, candidate_now.xulie(:, i-1))= Max;
%     candidate_now.xulie(i)= find(distemp(candidate_now.xulie(i-1), :)==min(distemp(candidate_now.xulie(i-1), :)), 1, 'first');
% end;
% figure(3); plot(X([Start,candidate_now.xulie,Start]), Y([Start,candidate_now.xulie,Start]), 'b-*');
% hold on
% plot(X(Start), Y(Start), 'r*'); hold off
% cat(1, candidate_now.xulie)
%====

candidate_now.value= dis(Start, candidate_now.xulie(1)); 
for i= 1:N-2
    candidate_now.value= candidate_now.value + dis(candidate_now.xulie(i), candidate_now.xulie(i+1));
end;
candidate_now.value= candidate_now.value + dis(candidate_now.xulie(end), Start); % 当前解（目标值）


% 初始化禁忌表
for i= 1:tabulength
    tabulist(i).change= [0, 0];
    tabulist(i).xulie= 0;
    tabulist(i).value= 0;
end;

% 目前最佳解
best_so_far.xulie= candidate_now.xulie;
best_so_far.value= candidate_now.value;
if displayflag
    figure(4); plot(X([Start,best_so_far.xulie,Start]), Y([Start,best_so_far.xulie,Start]), 'b-*'); 
end;

while runcount<T % 最多搜索100次
    runcount= runcount+1;
    TongJi1(runcount)= candidate_now.value; % 统计当前解的变化，用作显示
    TongJi2(runcount)= best_so_far.value; %
    %======== 生成邻域 获得候选集的城市序列和目标值
    kk= 1;
    for i= 1:N-2
        for j= i+1:N-1
            flag= true;
            for n= 1:tabulength
                tempflag= sum(abs(tabulist(n).change-[i, j]));
                flag= flag & tempflag;
            end;
            if flag
                temp= cat(1, candidate_now.xulie);
                %=== 两两城市互相交换顺序
                t= temp(i);
                temp(i)= temp(j);
                temp(j)= t;
                %===
                neighborhood(kk).xulie= temp; % 候选集的城市序列
                neighborhood(kk).change= [i, j];
                kk= kk+1;
            end;
        end;
    end;
    L= length(neighborhood);
    %=== 加入变异算子
    if IsBianYi && runcount<50
        i= L;
        for k= 1:ceil(N)
            if rand<Pby
                i= i+1;
                neighborhood(i).change= [0, 0];
                if rand<0.5
                    xulie_temp= best_so_far.xulie;
                else
                    xulie_temp= candidate_now.xulie;
                end;
                if rand<0.4 % 以一个城市为起始点，找到之后离他最近的一个城市，然后将他们之间的城市倒序
                    
                    if runcount<TT % TT之前用规则寻找起始点
                        disj= 0;
                        for j= 1:N-3
                            disj_temp= dis(xulie_temp(j), xulie_temp(j+1));
                            if disj_temp>disj
                                disj= disj_temp;
                                indexj= j;
                            end;
                        end;
                        if dis(xulie_temp(N-2), Start)<dis(xulie_temp(N-1), Start)
                            xulie_temp(end-1:end)= xulie_temp(end:-1:end-1);
                        end;
                    else % TT之后随机产生起始点
                        indexj= ceil(rand*(N-3));
                    end; % end if runcount<TT
%                     if indexj<N-3
                        disj= Max;
                        for j= indexj+2:N-1
                            disj_temp= dis(xulie_temp(indexj), xulie_temp(j));
                            if disj_temp<disj
                                disj= disj_temp;
                                indexjj= j;
                            end;
                        end;
                        xulie_temp(indexj+1:indexjj)= xulie_temp(indexjj:-1:indexj+1);
                        neighborhood(i).xulie= xulie_temp;
%                     end; % end if indexj<N-3
                elseif rand<0.6 % 以一个城市为分界点，将之前之后的城市交换
                    
                    indexjc= ceil(rand*(N-2));
                    temp11= xulie_temp(1:indexjc);
                    temp22= xulie_temp(indexjc+1:end);
                    xulie_temp= [temp22, temp11];
                    neighborhood(i).xulie= xulie_temp;
                else % 找一个起始点，一个长度，提取出这段城市；再在剩下的城市中找到离起始点城市最近的城市
                     % 然后将
                    
                    indexby2length= ceil(rand*ceil(N/4));
                    indexby2start1= ceil(rand*(N-1-indexby2length));
                    tempby21= xulie_temp(indexby2start1 : indexby2start1+indexby2length-1);
                    disby2= dis + Max*eye(N);
                    disby2(:, [1, tempby21])= Max;
                    minindexby2= find(disby2(xulie_temp(indexby2start1), :)==min(disby2(xulie_temp(indexby2start1), :)), 1, 'first');
                    xulie_temp(indexby2start1 : indexby2start1+indexby2length-1)= [];
                    indexby2start2= find(xulie_temp==minindexby2);
                    xulie_tempp= xulie_temp(1:indexby2start2);
                    if rand<0.5
                        xulie_tempp(indexby2start2+1:indexby2start2+1+indexby2length-1)= tempby21;
                    else
                        xulie_tempp(indexby2start2+1+indexby2length-1:-1:indexby2start2+1)= tempby21;
                    end;
                    xulie_tempp(indexby2start2+1+indexby2length:N-1)= xulie_temp(indexby2start2+1:end);
                    neighborhood(i).xulie= xulie_tempp;
                end; % end if rand<0.7
            end; % end if rand<Pby
        end; 
        L= i;
    end; % end if IsBianYi
    %===

    for i= 1:L
        d= dis(Start, neighborhood(i).xulie(1));
        for k= 1:N-2
            d= d + dis(neighborhood(i).xulie(k), neighborhood(i).xulie(k+1));
        end;
        d= d + dis(neighborhood(i).xulie(end), Start);
        neighborhood(i).value= d; % 候选集的目标值
    end;
    %========
    candidate_next= sort(cat(1, neighborhood.value), 1).'; % 将后选集的目标值按从小到大排序
    if candidate_next(1)<best_so_far.value % 满足藐视准则
        index= find(cat(1, neighborhood.value)==candidate_next(1)); % 找到最小目标值对应的对象的位置
        L= length(index); % 有可能不只一个最小值
        if L>1 % 如果不只一个最小值 选择后面一种
            candidate_now.xulie= neighborhood(index(L)).xulie;
            candidate_now.change= neighborhood(index(L)).change;
        else % 只有一个最小值
            candidate_now.xulie= neighborhood(index(1)).xulie;
            candidate_now.change= neighborhood(index(L)).change;
        end;
        candidate_now.value= candidate_next(1);
        %=== 将这个对象放入禁忌表
        for i= tabulength-1:-1:1 % i= 2, 1
            tabulist(i+1).change= tabulist(i).change;
            tabulist(i+1).xulie= tabulist(i).xulie;
            tabulist(i+1).value= tabulist(i).value;
        end;
        tabulist(1).change= candidate_now.change;
        tabulist(1).xulie= candidate_now.xulie;
        tabulist(1).value= candidate_now.value;
        %===
        % 目前最佳解 更新
        best_so_far.xulie= candidate_now.xulie;
        best_so_far.value= candidate_now.value;
        if displayflag
            figure(4); plot(X([Start,best_so_far.xulie,Start]), Y([Start,best_so_far.xulie,Start]), 'b-*');
            hold on
            plot(X(Start), Y(Start), 'r*');
            hold off
        end;
    else % 不满足藐视规则，则选取最小的前5个中的非禁忌对象的最小一个作为当前解
        figure(4);
        for n= 1:10
            index= find(cat(1, neighborhood.value)==candidate_next(n));
            tx= neighborhood(index(1)).xulie;
            tc= neighborhood(index(1)).change;
            %=== 判断是否在禁忌表中
            flag= true;
            for i= 2:tabulength
                flag= flag & sum(abs(tc - tabulist(i).change))>0;
            end;
            %====
            if flag
                candidate_now.xulie= tx;
                candidate_now.change= tc;
                candidate_now.value= candidate_next(n);
                %=== 将这个对象放入禁忌表
                for i= tabulength-1:-1:1 % i= 2, 1
                    tabulist(i+1).change= tabulist(i).change;
                    tabulist(i+1).xulie= tabulist(i).xulie;
                    tabulist(i+1).value= tabulist(i).value;
                end;
                tabulist(1).change= candidate_now.change;
                tabulist(1).xulie= candidate_now.xulie;
                tabulist(1).value= candidate_now.value;
                %===
                break;
            end;
        end;
    end;

    clear neighborhood
    disp([num2str(runcount), '  ', num2str(cat(1, candidate_now.xulie), '%.0f'), '  ', num2str(cat(1, candidate_now.value), '%.5f'), '  ', num2str(cat(1, best_so_far.value), '%.5f')]); % 显示中间序列结果
end;

result= best_so_far.value; % 返回目标值结果
shunxu= best_so_far.xulie; % 返回城市顺序结果

%=== 显示结果
if displayflag
    figure(4), close;
    figure(4); plot(X([Start,shunxu,Start]), Y([Start,shunxu,Start]), 'b*'); 
    figure(4), hold on
    plot(X(Start), Y(Start), 'r*');
    plot([X(Start),X(shunxu(1))], [Y(Start),Y(shunxu(1))]);
    pause(0.5)
    for k= 1:N-2
        figure(4), plot([X(shunxu(k)),X(shunxu(k+1))], [Y(shunxu(k)),Y(shunxu(k+1))]);
        pause(0.5)
    end;
    figure(4),plot([X(shunxu(end)),X(Start)], [Y(shunxu(end)),Y(Start)]);
    pause(0.5)
    figure(4), hold off
    
    figure(2), plot(1:length(TongJi1), TongJi1, 'b-', 1:length(TongJi2), TongJi2, 'r-'); legend('now', 'best so far'); title('process of refreshing solution'); xlabel('times'); ylabel('distance')

    disp('结果：')
    disp(['迭代次数：', num2str(runcount)]);
    disp(['最短路程：', num2str(best_so_far.value)]);
    disp(['路线：', num2str([Start, best_so_far.xulie, Start])]);

    disp('验证：')
    d= dis(Start,best_so_far.xulie(1));
    for k= 1:N-2
        d= d + dis(best_so_far.xulie(k), best_so_far.xulie(k+1));
    end;
    d= d + dis(best_so_far.xulie(end), Start);

    disp(['验证路程：', num2str(d)]);

    disp('@CopyRight Cheng Defu')
end;
%===

return;