禁忌搜索算法代码注释

在网上看到了一位别人写的禁忌搜索的代码,看了半天才理解,为了让大家更好地理解,我在别人代码的基础上添加了详细的注释,代码我按照算法流程图注释的,配合流程图看会更容易理解。大家可以把代码复制到matlab里看,有详细的分块,看起来会清楚一些

这里是原文:禁忌搜索算法(Tabu Search)

在这里插入图片描述

%% 一个旅行商人要拜访全国31个省会城市,且每个城市只能拜访一次,求所有路径之中的最小值
% 禁忌搜索算法求解TSP问题
function [BestShortcut,theMinDistance]=TabuSearch
clear;
clc;
% 全国31个省会城市坐标
Clist=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;3238 1229;...
    4196 1044;4312  790;4386  570;3007 1970;2562 1756;2788 1491;2381 1676;...
    1332  695;3715 1678;3918 2179;4061 2370;3780 2212;3676 2578;4029 2838;...
    4263 2931;3429 1908;3507 2376;3394 2643;3439 3201;2935 3240;3140 3550;...
    2545 2357;2778 2826;2370 2975];

%% 计算距离矩阵
CityNum=size(Clist,1); % TSP问题的规模,即城市数目
dislist=zeros(CityNum);
for i=1:CityNum
    for j=1:CityNum
        dislist(i,j)=((Clist(i,1)-Clist(j,1))^2+(Clist(i,2)-Clist(j,2))^2)^0.5;
    end
end
%% 定义属性
TabuList=zeros(CityNum);                      % 禁忌表设置为所有城市互换对,初始禁忌长度都为0,即禁忌表为空
TabuLength=round((CityNum*(CityNum-1)/2)^0.5);% 禁忌表长度(tabu length),设置为二领域个数的开方,round四舍五入
CandidatesNum=200;                            % 领域集的个数 (全部二领域解个数太多,我们随机选取200个),由这些领域解中选出候选解
Candidates=zeros(CandidatesNum,CityNum);      % 选取的领域解集合,每一行代表一个领域解,由这些领域解中选出候选解
S0=randperm(CityNum);                         % 随机产生初始解,返回一行1到CityNum的整数数组
BSF=S0;                                       % 当前最佳路径
BestL=Inf;                                    % 当前最佳解距离,inf为极大值
p=1;                                          % 记录迭代次数
StopL=1000;                                   % 最大迭代次数
%% 用户控制界面
figure(1);
stop = uicontrol('style','toggle','string'...
    ,'stop','background','white');
tic;                                          % 用来保存当前时间
%% 禁忌搜索循环
while p<StopL
    if CandidatesNum>CityNum*(CityNum-1)/2
        disp('领域解个数不大于n*(n-1)/2!');
        break;
    end
    
    ALong(p)=DistanceOfS(dislist,S0);       % 当前解适配值,也就是当前路径的长度
    
    %% 以下while循环生成随机的[200,2]的矩阵A,每一个元素都是在1-31之间的
    % 用于后面为当前解生成邻域解
    i=1;
    A=zeros(CandidatesNum,2);       % 解中交换的城市矩阵
    while i<=CandidatesNum
        M=CityNum*rand(1,2); 
        M=ceil(M);                  % 返回一行两列数值为 1-CityNum 之间的数组
        %【1-----------------------------
        if M(1)~=M(2)
            A(i,1)=max(M(1),M(2));
            A(i,2)=min(M(1),M(2));
            %【2-------------------------
            if i==1
                isa=0;
            else
                %【3---------------------
                % 判断是否生成相同的交换城市组,如果是,重新生成
                for j=1:i-1
                    if A(i,1)==A(j,1) && A(i,2)==A(j,2)
                        isa=1;
                        break;
                    else
                        isa=0;
                    end
                end
                %3】---------------------
            end
            %2】-------------------------
            %【4---------
            if ~isa
                i=i+1;
            end
            %4】---------
        end
        %1】-----------------------------
    end
    %% 产生领域解,选取前100个为确定候选解
    % 保留前BestCandidateNum个最好领域解作为候选解
    BestCandidateNum=100;
    % 四列,分别存数i、路径长度、交换城市A(i,1)、交换城市A(i,2)
    BestCandidate=Inf*ones(BestCandidateNum,4); 
    F=zeros(1,CandidatesNum);
    
    % 这相当于是产生一个S0的邻域...包含CandidatesNum个领域解,最后选出的BestCandidate才是候选解
    for i=1:CandidatesNum
        Candidates(i,:)=S0;  % 领域解集合
        Candidates(i,[A(i,2),A(i,1)])=S0([A(i,1),A(i,2)]);
        % 计算当前解适配值,也就是当前路径的长度
        F(i)=DistanceOfS(dislist,Candidates(i,:));
        %----------------------------------
        % 选出100个最好的领域解
        if i<=BestCandidateNum
            BestCandidate(i,2)=F(i);
            BestCandidate(i,1)=i;
            BestCandidate(i,3)=S0(A(i,1));
            BestCandidate(i,4)=S0(A(i,2));
        else
            for j=1:BestCandidateNum
                if F(i)<BestCandidate(j,2)
                    BestCandidate(j,2)=F(i);
                    BestCandidate(j,1)=i;
                    BestCandidate(j,3)=S0(A(i,1));
                    BestCandidate(j,4)=S0(A(i,2));
                    break;
                end
            end
        end
        %----------------------------------
    end
    % 对BestCandidate排序
    [JL,Index]=sort(BestCandidate(:,2)); % JL是排序后的向量(用不到),index是JL中每一项对应于BestCandidate(:,2)中项的索引,排序是按升序进行的。
    SBest=BestCandidate(Index,:);
    BestCandidate=SBest;                 % 此时的BestCandidate是按第二列路径长度升序排列的
    %% 特赦准则是否满足,更新禁忌表
    % 此时BestCandidate(1,1)存储的i下标代表Candidates中最短路径对应的行下标
    %% Y
    if BestCandidate(1,2)<BestL
        BestL=BestCandidate(1,2);
        % 当前解(路经集合)变成候选集合的第BestCandidate(1,1)行,这一行最小
        S0=Candidates(BestCandidate(1,1),:);
        BSF=S0;  
        for m=1:CityNum
            for n=1:CityNum
                if TabuList(m,n)~=0
                    TabuList(m,n)=TabuList(m,n)-1;                         % 更新禁忌表,禁忌长度减1
                end
            end
        end
        TabuList(BestCandidate(1,3),BestCandidate(1,4))=TabuLength;        % 更新禁忌表,设置当前解的禁忌长度
    %% N
    % 此时没有候选解比当前解更优
    else 
        for i=1:BestCandidateNum
            if  TabuList(BestCandidate(i,3),BestCandidate(i,4))==0         % 如果已经解禁
                S0=Candidates(BestCandidate(i,1),:);
                for m=1:CityNum
                    for n=1:CityNum
                        if TabuList(m,n)~=0
                            TabuList(m,n)=TabuList(m,n)-1;                 % 更新禁忌表,禁忌长度减1
                        end
                    end
                end
                TabuList(BestCandidate(i,3),BestCandidate(i,4))=TabuLength;% 更新禁忌表,设置当前解的禁忌长度
                break;                                                     % 找到第一个非禁忌的最好的候选解作为当前解,退出for循环
            end
        end
    end
    % 更新当前最好的解适配值,也就是当前最短路径的长度
    % BestCandidate(1,2)<BestL时更新,否则不更新,放在这里是为了方便后面画适应度曲线用
    ArrBestL(p)=BestL;
    %% 每次迭代都画出最好的路径
    for i=1:CityNum-1
        plot([Clist(BSF(i),1),Clist(BSF(i+1),1)],[Clist(BSF(i),2),Clist(BSF(i+1),2)],'bo-');    % 画出两点的连线
        hold on;
    end
    plot([Clist(BSF(CityNum),1),Clist(BSF(1),1)],[Clist(BSF(CityNum),2),Clist(BSF(1),2)],'ro-');% 画出最后一点与出发点的连线(红色)
    title(['迭代次数:',int2str(p),' 优化最短距离:',num2str(BestL)]);
    hold off;
    pause(0.005);
    if get(stop,'value')==1
        break;
    end
    % 存储中间结果为图片
    if (p==1||p==5||p==10||p==20||p==60||p==150||p==400||p==800||p==1500||p==2000)
        filename=num2str(p);
        fileformat='jpg';
        saveas(gcf,filename,fileformat);
    end
    p=p+1;                                                                 % 迭代次数加1
end
%% 
toc;                                         % 用来保存完成时间
BestShortcut=BSF;                            % 最佳路线
theMinDistance=BestL;                        % 最佳路线长度
set(stop,'style','pushbutton','string','close', 'callback','close(gcf)');
%% 画适应度进化曲线
figure(2);
plot(ArrBestL,'b');
xlabel('迭代次数');
ylabel('目标函数值');
title('适应度进化曲线');
grid;
hold on;
%% 运行时间
% figure(3)
% plot(toc-tic,'b');
% grid;
% title('运行时间');
% legend('Best So Far','当前解');
% 计算适配值函数,即路径的长度
function F=DistanceOfS(dislist,s)   
DistanV=0;
n=size(s,2);
for i=1:(n-1)
    DistanV=DistanV+dislist(s(i),s(i+1));
end
    DistanV=DistanV+dislist(s(n),s(1));      
F=DistanV;

运行结果:
在这里插入图片描述
在这里插入图片描述

  • 28
    点赞
  • 107
    收藏
    觉得还不错? 一键收藏
  • 4
    评论
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值