蚁群算法-求解TSP问题

蚁群算法(ant colony optimization, ACO)

1 原理介绍

​ 蚁群算法（Ant Colony Algorithm）最初于1992年由意大利学者M.Dorigo等人提出，它是一种模拟自然界中真实蚁群觅食行为的仿生优化算法。研究发现：每只蚂蚁觅食时在走过的路线上会留下一种称为信息素的物质，蚂蚁之间靠感知这种物质的浓度进行信息传递。蚂蚁在选择路径时总是倾向于朝信息索浓度高的方向移动，而距离短的路径上走过的蚂蚁多，留下的信息素也多，后续蚂蚁选择它的概率也会越大；其他路径上的信息素会随着时间的推移不断挥发，这样就形成了一种正反馈机制，最后整个蚁群聚集到最短路径上。

​ 人工蚁群算法模拟了这一过程。每只蚂蚁在解空间独立地搜索可行解，解越好留下的信息素越多，随着算法推进，较优解路径上的信息素增多，选择它的蚂蚁也随之增多，最终收敛到最优或近似最优的解上。

 

 

 

2 代码实现

2.1 MATALB实现：

%%%一个旅行商人要拜访全国31个省会城市，需要选择最短的路径%%%%

%%%蚁群算法解决TSP问题%%%%%%%
 
clear all; %清除所有变量
close all; %清图
clc;       %清屏

m=50;       %% m 蚂蚁个数
Alpha=1;    %% Alpha 表征信息素重要程度的参数
Beta=5;     %% Beta 表征启发式因子重要程度的参数
Rho=0.1;    %% Rho 信息素蒸发系数
NC_max=100; %%最大迭代次数
Q=100;      %%信息素增加强度系数

C=[
1304 2312;
3639 1315;
4177 2244;
3712 1399;
3488 1535;
3326 1556;
3238 1229;
4196 1004;
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 2367;
3394 2643;
3439 3201;
2935 3240;
3140 3550;
2545 2357;
2778 2826;
2370 2975
];%%31个省会坐标
%%-------------------------------------------------------------------------
%% 主要符号说明
%% C n个城市的坐标，n×2的矩阵
%% NC_max 最大迭代次数
%% m 蚂蚁个数
%% Alpha 表征信息素重要程度的参数
%% Beta 表征启发式因子重要程度的参数
%% Rho 信息素蒸发系数
%% Q 信息素增加强度系数
%% R_best 各代最佳路线
%% L_best 各代最佳路线的长度
%%=========================================================================
%% 第一步：变量初始化
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; %i=j时不计算，应该为0，但后面的启发因子要取倒数，用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(NC_max,n);       %各代最佳路线
L_best=inf.*ones(NC_max,1);   %各代最佳路线的长度
L_ave=zeros(NC_max,1);        %各代路线的平均长度

while NC<=NC_max  %停止条件之一：达到最大迭代次数，停止
    %%第二步：将m只蚂蚁放到n个城市上
    Randpos=[];   %随即存取
    for i=1:(ceil(m/n))
        Randpos=[Randpos,randperm(n)]; %执行两次操作之后Randpos的维度为1*62
    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                            %找到未访问的城市，并存在数组J中
                if isempty(find(visited==k, 1))  %开始时置0 find函数返回在visited数组中k所在的位置 没有则返回0 1表示只找1次
                    J(Jc)=k;
                    Jc=Jc+1;                     %访问的城市个数自加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);         %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,:);%保留一下上次最优路线至tabu第一行，保障本次迭代情况不至于太差
    end
    
    %%第四步：记录本次迭代每只蚂蚁所走距离L，记录每次迭代最佳路线距离L_best和最佳路线信息R_best
    L=zeros(m,1);                       %开始距离为0，m*1的列向量
    for i=1:m
        R=Tabu(i,:);
        for j=1:(n-1)
            L(i)=L(i)+D(R(j),R(j+1));   %原距离加上第j个城市到第j+1个城市的距离
        end
        L(i)=L(i)+D(R(1),R(n));         %一轮下来后走过的距离
    end
    
    L_best(NC)=min(L);                  %最佳距离取最小
    L_ave(NC)=mean(L);                  %此轮迭代后的平均距离
    
    pos=find(L==L_best(NC));
    R_best(NC,:)=Tabu(pos(1),:);        %此轮迭代后的最佳路线
    
    NC=NC+1                             %迭代继续

    %%第五步：更新信息素
    Delta_Tau=zeros(n,n);               %开始时信息素为n*n的0矩阵
    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);
            %此次循环在路径（i，j）上的信息素增量
        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);            %直到最大迭代次数
end

%%第七步：输出结果
Pos=find(L_best==min(L_best)); %找到最佳路径（非0为真）
Shortest_Route=R_best(Pos(1),:)%最大迭代次数后最佳路径
Shortest_Length=L_best(Pos(1)) %最大迭代次数后最短距离

%% 画出路线图，和L_best，L_ave迭代曲线
figure(1)
subplot(1,2,1)                  %绘制第一个子图形
N=length(Shortest_Route);
scatter(C(:,1),C(:,2));
for i = 1:size(C,1)
	text(C(i,1),C(i,2),['   ' num2str(i)]);
end
hold on
plot([C(Shortest_Route(1),1),C(Shortest_Route(N),1)],[C(Shortest_Route(1),2),C(Shortest_Route(N),2)],'g')
hold on
for ii=2:N
    plot([C(Shortest_Route(ii-1),1),C(Shortest_Route(ii),1)],[C(Shortest_Route(ii-1),2),C(Shortest_Route(ii),2)],'g')
    hold on
end
title('旅行商问题优化结果 ')

subplot(1,2,2)                  %绘制第二个子图形
plot(L_best)
hold on                         %保持图形
plot(L_ave,'r')
title('平均距离和最短距离')      %标题

在MATLAB的算法实现当中，定义蚂蚁个数m=50，城市个数n=31，信息素重要程度参数Alpha=1，启发因子重要程度Beta=31，信息素蒸发系数Rho=0.1，最大迭代次数NC_max=200，信息素增强系数Q=100。定义禁忌表Tabu(50,31)，启发因子表Eta(31,31)，信息素浓度表Tau(31,31)，蚂蚁选择下一个未走城市概率P(k)。

启发因子表Eta(31,31)和信息素浓度表Tau(31,31)共同决定蚂蚁选择下一个未走城市概率P(k)，其计算过程为（2）式。在这之中，禁忌表Tabu(50,31)在每只蚂蚁每走过一个城市都进行更新，启发因子表Eta(31,31)为城市之间的距离始终保持不变，信息素浓度表Tau(31,31)在禁忌表填满一次后根据（1）式进行更新。

2.2 C++实现：

2.3 Python实现：

import numpy as np
import matplotlib.pyplot as plt
import pylab

coordinates = np.array([[0.8223865, 0.90249145],
                        [0.87299287, 0.97785658],
                        [0.58388132, 0.31408447],
                        [0.72751158, 0.05415505],
                        [0.60553193, 0.00697702],
                        [0.45564878, 0.15191931],
                        [0.411461, 0.17028803],
                        [0.42169505, 0.29723746],
                        [0.34342426, 0.13354594],
                        [0.26429325, 0.17590559],
                        [0.19763023, 0.12130362],
                        [0.13727432, 0.25075126],
                        [0.00319099, 0.38633692],
                        [0.24217939, 0.40168494],
                        [0.21614921, 0.42319542],
                        [0.26362807, 0.54363507],
                        [0.0312916, 0.64368314],
                        [0.36832452, 0.75778461],
                        [0.50299907, 0.50413817],
                        [0.65584064, 0.93553257]])


def getdistmat(coordinates):
    num = coordinates.shape[0]
    distmat = np.zeros((20, 20))
    for i in range(num):
        for j in range(i, num):
            distmat[i][j] = distmat[j][i] = np.linalg.norm(coordinates[i] - coordinates[j])
    return distmat


distmat = getdistmat(coordinates)
numant = 40  # 蚂蚁个数
numcity = coordinates.shape[0]  # 城市个数
alpha = 1  # 信息素重要程度因子
beta = 5  # 启发函数重要程度因子
rho = 0.1  # 信息素的挥发速度
Q = 1
iter = 0
itermax = 250
etatable = 1.0 / (distmat + np.diag([1e10] * numcity))  # 启发函数矩阵，表示蚂蚁从城市i转移到矩阵j的期望程度
pheromonetable = np.ones((numcity, numcity))  # 信息素矩阵
pathtable = np.zeros((numant, numcity)).astype(int)  # 路径记录表
distmat = getdistmat(coordinates)  # 城市的距离矩阵
lengthaver = np.zeros(itermax)  # 各代路径的平均长度
lengthbest = np.zeros(itermax)  # 各代及其之前遇到的最佳路径长度
pathbest = np.zeros((itermax, numcity))  # 各代及其之前遇到的最佳路径长度

while iter < itermax:
    # 随机产生各个蚂蚁的起点城市
    if numant <= numcity:  # 城市数比蚂蚁数多
        pathtable[:, 0] = np.random.permutation(range(0, numcity))[:numant]
    else:  # 蚂蚁数比城市数多，需要补足
        pathtable[:numcity, 0] = np.random.permutation(range(0, numcity))[:]
        pathtable[numcity:, 0] = np.random.permutation(range(0, numcity))[:numant - numcity]
    length = np.zeros(numant)  # 计算各个蚂蚁的路径距离
    for i in range(numant):
        visiting = pathtable[i, 0]  # 当前所在的城市
        unvisited = set(range(numcity))  # 未访问的城市,以集合的形式存储{}
        unvisited.remove(visiting)  # 删除元素；利用集合的remove方法删除存储的数据内容
        for j in range(1, numcity):  # 循环numcity-1次，访问剩余的numcity-1个城市
            # 每次用轮盘法选择下一个要访问的城市
            listunvisited = list(unvisited)
            probtrans = np.zeros(len(listunvisited))
            for k in range(len(listunvisited)):
                probtrans[k] = np.power(pheromonetable[visiting][listunvisited[k]], alpha) \
                               * np.power(etatable[visiting][listunvisited[k]], alpha)
            cumsumprobtrans = (probtrans / sum(probtrans)).cumsum()
            cumsumprobtrans -= np.random.rand()
            k = listunvisited[(np.where(cumsumprobtrans > 0)[0])[0]]  # python3中原代码运行bug，类型问题；鉴于此特找到其他方法
            # 通过where（）方法寻找矩阵大于0的元素的索引并返回ndarray类型，然后接着载使用[0]提取其中的元素，用作listunvisited列表中
            # 元素的提取（也就是下一轮选的城市）
            pathtable[i, j] = k  # 添加到路径表中（也就是蚂蚁走过的路径)
            unvisited.remove(k)  # 然后在为访问城市set中remove（）删除掉该城市
            length[i] += distmat[visiting][k]
            visiting = k
        length[i] += distmat[visiting][pathtable[i, 0]]  # 蚂蚁的路径距离包括最后一个城市和第一个城市的距离
        # 包含所有蚂蚁的一个迭代结束后，统计本次迭代的若干统计参数
    lengthaver[iter] = length.mean()
    if iter == 0:
        lengthbest[iter] = length.min()
        pathbest[iter] = pathtable[length.argmin()].copy()
    else:
        if length.min() > lengthbest[iter - 1]:
            lengthbest[iter] = lengthbest[iter - 1]
            pathbest[iter] = pathbest[iter - 1].copy()
        else:
            lengthbest[iter] = length.min()
            pathbest[iter] = pathtable[length.argmin()].copy()
    # 更新信息素
    changepheromonetable = np.zeros((numcity, numcity))
    for i in range(numant):
        for j in range(numcity - 1):
            changepheromonetable[pathtable[i, j]][pathtable[i, j + 1]] += Q / distmat[pathtable[i, j]][
                pathtable[i, j + 1]]  # 计算信息素增量
        changepheromonetable[pathtable[i, j + 1]][pathtable[i, 0]] += Q / distmat[pathtable[i, j + 1]][pathtable[i, 0]]
    pheromonetable = (1 - rho) * pheromonetable + changepheromonetable  # 计算信息素公式
    iter += 1  # 迭代次数指示器+1
    print("iter:", iter)

# 做出平均路径长度和最优路径长度
fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(12, 10))
axes[0].plot(lengthaver, 'k', marker=u'')
axes[0].set_title('Average Length')
axes[0].set_xlabel(u'iteration')

axes[1].plot(lengthbest, 'k', marker=u'')
axes[1].set_title('Best Length')
axes[1].set_xlabel(u'iteration')
fig.savefig('average_best.png', dpi=500, bbox_inches='tight')
plt.show()
print(lengthbest[-1])

# 作出找到的最优路径图
bestpath = pathbest[-1]
plt.plot(coordinates[:, 0], coordinates[:, 1], 'r.', marker=u'$\cdot$')
plt.xlim([0, 1])
plt.ylim([0, 1])

for i in range(numcity - 1):
    m = int(bestpath[i])
    n = int(bestpath[i + 1])
    plt.plot([coordinates[m][0], coordinates[n][0]], [coordinates[m][1], coordinates[n][1]], 'k')
plt.plot([coordinates[int(bestpath[0])][0], coordinates[int(n)][0]],
         [coordinates[int(bestpath[0])][1], coordinates[int(n)][1]], 'b')
ax = plt.gca()
ax.set_title("Best Path, length is % s" % lengthbest[-1])
ax.set_xlabel('X axis')
ax.set_ylabel('Y_axis')

plt.savefig('best path.png', dpi=500, bbox_inches='tight')
plt.show()