从Matlab实例学习蚁群算法

前言

个人学习完蚁群算法，觉得可以用一句话来概括其本质核心：

挑选可行解进行遍历， 再根据结果更新下一轮可行解的挑选概率。

有点类似于强化学习， 更好的游戏得分会进一步鼓励下一次类似的尝试。 只是换了一种生物学的说法，这种鼓励对应于蚂蚁分泌的信息素——信息素越多的路径会吸引更多蚂蚁探索。 目前来看，所有的各类仿生启发式算法，本质上仍始终是试图达到比穷搜更高效的搜索， 在缺乏数学理论解释的情况下， 便以生物学作为启发的基石。

蚁群算法求解TSP问题

TSP指的是旅行商问题， 简单而言， 有$N$个城市， 需要规划一条路径，经过所有城市均一次， 且总路径最短。 蚁群算法的步骤如下：

1. 将 $m$ 只蚂蚁各自随机放到一个城市之中。 因为我们要试图构造$m$个可行解，因此初始点先随意选在一个城市上。这一步是纯随机选取。
2. $m$ 只蚂蚁各自按照概率选择下一个城市， 直到走完所以所有城市， 每个城市只能走一次。 每次出发自城市$i$， 下一步选择城市$j$的概率为：
   $$p_{ij}\sim {\tau }_{ij}{\eta }_{ij}\text{\hspace{1em}(1)}$$
   其中， $\tau$代表信息素因子，初始化全为$1$， 正比于尝试过这一路径的蚂蚁数； $\eta$代表启发式因子，在本问题中，反比于城市$i$与城市$j$之间的距离， 因为我们问题的优化目标是最小化总路程。
   $m$只蚂蚁均完成周游后， 我们相当于得到了$m$个TSP问题的可行解。
3. 计算这$m$个可行解的结果，并记录最佳路径。
4. 根据本次结果，更新信息素信息：
   $${\tau }_{ij}=(1-\rho ){\tau }_{ij}+\Delta {\tau }_{ij}\text{\hspace{1em}(2)}$$
   其中$\rho$代表蒸发系数， 对应于自然界中，蚂蚁留下信息素的蒸发。 $\Delta {\tau }_{ij}$则代表了基于此次蚂蚁探索结果的更新量， 正比于本次经过了从城市$i$到城市$j$这一路径的蚂蚁数量， 具体的关系则不同算法准则间有所区别。 但本质都是，被更多蚂蚁探索过的解空间，会进一步鼓励以后的蚂蚁继续探索。
5. 蚂蚁搜寻代数是否达到预设值。 若否，跳转到步骤1继续探索。记录每次探索的最佳路径，最终输出所有探索的最佳结果。

简单总结： 基于过往探索经验 + 基于目标的启发式因子， 优化未来的搜索。

一些改进

基于蚁群算法的特点， 很容易想到一些改进方案

• 精英蚂蚁： 对应于更好结果的蚂蚁，释放更多的信息素——也即该蚂蚁的信息素权重更高。
• 最大最小蚂蚁： 1. 每次循环后只有一只蚂蚁能更新信息素 2. 每个解元素(如TSP问题中的每两个城市的转移)上的信息素有上下界限制。 3. 算法初始阶段将信息素重新初始化。 目的：防止某个解元素的信息素不断增加，阻碍了对其他解的探索。
• 排序蚂蚁： 精英蚂蚁的细化升级版。 每个蚂蚁分泌信息素的权重正比于其寻找到的解的性能好坏。
• 自适应蚁群： 主要两个改动 1. 每次找出最优解并保留。 2. 自适应改变蒸发系数$\rho$， 防止由于蒸发过快使得一些从未被搜索过的路径信息素越来越少，再也不被探测。 因此设定当算法在多代搜索后仍没有提升时， 降低蒸发系数，鼓励新的探索。

蚁群算法解决TSP问题Matlab代码

%%%%%%%%%%%%%%%%%%%%蚁群算法解决TSP问题%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;                %清除所有变量
close all;                %清图
clc;                      %清屏
m=50;                     %蚂蚁个数
Alpha=1;                  %信息素重要程度参数              
Beta=5;                   %启发式因子重要程度参数
Rho=0.1;                  %信息素蒸发系数
G_max=200;                %最大迭代次数
Q=100;                    %信息素增加强度系数
C=[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];                 %31个省会城市坐标
%%%%%%%%%%%%%%%%%%%%%%%%第一步：变量初始化%%%%%%%%%%%%%%%%%%%%%%%%
n=size(C,1);              %n表示问题的规模（城市个数）
D=zeros(n,n);             %D表示两个城市距离间隔矩阵
for i=1:n
    for j=1:n
        if i~=j
            D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;
        else
            D(i,j)=eps;
        end
        D(j,i)=D(i,j);
    end
end
Eta=1./D;                    %Eta为启发因子，这里设为距离的倒数
Tau=ones(n,n);               %Tau为信息素矩阵
Tabu=zeros(m,n);             %存储并记录路径的生成
NC=1;                        %迭代计数器
R_best=zeros(G_max,n);       %各代最佳路线
L_best=inf.*ones(G_max,1);   %各代最佳路线的长度
figure(1);%优化解
while NC<=G_max            
    %%%%%%%%%%%%%%%%%%第二步：将m只蚂蚁放到n个城市上%%%%%%%%%%%%%%%%
    Randpos=[];
    for i=1:(ceil(m/n))
        Randpos=[Randpos,randperm(n)];
    end
    Tabu(:,1)=(Randpos(1,1:m))'; 
    %%%%%第三步：m只蚂蚁按概率函数选择下一座城市，完成各自的周游%%%%%%
    for j=2:n
        for i=1:m
            visited=Tabu(i,1:(j-1));  %已访问的城市
            J=zeros(1,(n-j+1));       %待访问的城市
            P=J;                      %待访问城市的选择概率分布
            Jc=1;
            for k=1:n
                if length(find(visited==k))==0
                    J(Jc)=k;
                    Jc=Jc+1;
                end
            end
            %%%%%%%%%%%%%%%%%%计算待选城市的概率分布%%%%%%%%%%%%%%%%
            for k=1:length(J)
                P(k)=(Tau(visited(end),J(k))^Alpha)...
                    *(Eta(visited(end),J(k))^Beta);
            end
            P=P/(sum(P));
            %%%%%%%%%%%%%%%%按概率原则选取下一个城市%%%%%%%%%%%%%%%%
            Pcum=cumsum(P);
            Select=find(Pcum>=rand);
            to_visit=J(Select(1));
            Tabu(i,j)=to_visit;
        end
    end
    if NC>=2
        Tabu(1,:)=R_best(NC-1,:);
    end
    %%%%%%%%%%%%%%%%%%%第四步：记录本次迭代最佳路线%%%%%%%%%%%%%%%%%%
    L=zeros(m,1);
    for i=1:m
        R=Tabu(i,:);
        for j=1:(n-1)
            L(i)=L(i)+D(R(j),R(j+1));
        end
        L(i)=L(i)+D(R(1),R(n));
    end
    L_best(NC)=min(L);
    pos=find(L==L_best(NC));
    R_best(NC,:)=Tabu(pos(1),:);
    %%%%%%%%%%%%%%%%%%%%%%%%%第五步：更新信息素%%%%%%%%%%%%%%%%%%%%%%
    Delta_Tau=zeros(n,n);
    for i=1:m
        for j=1:(n-1)
            Delta_Tau(Tabu(i,j),Tabu(i,j+1))=...
                Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i);
        end
        Delta_Tau(Tabu(i,n),Tabu(i,1))=...
            Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i);
    end
    Tau=(1-Rho).*Tau+Delta_Tau;
    %%%%%%%%%%%%%%%%%%%%%%%第六步：禁忌表清零%%%%%%%%%%%%%%%%%%%%%%
    Tabu=zeros(m,n);
    %%%%%%%%%%%%%%%%%%%%%%%%%历代最优路线%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:n-1
        plot([ C(R_best(NC,i),1), C(R_best(NC,i+1),1)],...
            [C(R_best(NC,i),2), C(R_best(NC,i+1),2)],'bo-');
        hold on;
    end
    plot([C(R_best(NC,n),1), C(R_best(NC,1),1)],...
        [C(R_best(NC,n),2), C(R_best(NC,1),2)],'ro-');  
    title(['优化最短距离:',num2str(L_best(NC))]);
    hold off;
    pause(0.005);
    NC=NC+1;    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%第七步：输出结果%%%%%%%%%%%%%%%%%%%%%%%%%%
Pos=find(L_best==min(L_best));
Shortest_Route=R_best(Pos(1),:);            %最佳路线
Shortest_Length=L_best(Pos(1));             %最佳路线长度
figure(2),
plot(L_best)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120