39行代码实现蚁群算法解决TSP问题程序（python）

见过不少用蚁群算法解决TSP问题的程序，多数都要过百行代码，显得冗长，而且效率低。
我这里也提供一个短小的程序。全部程序，除去空行和注释行，仅39行代码。
同时采用普遍撒微量信息素的方式，避免发生彻底的收敛，保证蚂蚁有概率突破局部最优解。同时也避免了分母过小，被计算机认为是零而报错。
本程序可获得非常接近最优的解：对dantzig42问题可以得到最优解：699；对环游中国34城市问题大概率可以获得16000公里以内的解。

import numpy as np

#准备工作
#读取城市矩阵
distance_matrix=np.load(r'distance_matrix.npy')#提前准备的距离矩阵文件
global city_num,ants_number,best_route,best,tao#城市数,蚁数，已知最优解，最优距离，信息素矩阵，为全局变量
city_num=len(distance_matrix)
ants_number=city_num

#距离矩阵变型
d0=1000/distance_matrix#距离矩阵的倒数*1000
for i in range(city_num):d0[i,i]=0

#计算一个解的总路程函数
def distance_total(route):
    return distance_matrix[route[np.arange(-1,city_num-1)],route].sum()

#初始最佳路径
best_route=np.random.permutation(city_num)
best=distance_total(best_route)

#初始信息素矩阵
tao=1-np.identity(city_num)#全为1

#一只蚂蚁走一步
def onestep(trace):
    to_be_chosen=np.setdiff1d(np.arange(city_num),trace)
    d1=(tao*(d0**2))[trace[-1],to_be_chosen]
    p=d1/sum(d1)#概率
    return np.random.choice(to_be_chosen,size=1,p=p)

result_this=[];result_best=[]#绘图用
Generations=3000#迭代3000代后终止
for gen in range(Generations):
    #所有蚂蚁走一圈回到出发点
    ants_matrix=np.random.randint(0,city_num-1,size=(ants_number,1))#初始随机布置蚂蚁
    for i in range(city_num-1):
        ants_matrix=np.concatenate((ants_matrix,np.apply_along_axis(onestep,arr=ants_matrix,axis=1)),axis=1)#所有走一圈
    distance_all=np.apply_along_axis(distance_total,arr=ants_matrix,axis=1)#计算每只蚂蚁走过的距离
    if distance_all.min()<best:
        best_route=ants_matrix[distance_all.argmin()]#最短路径
        best=distance_all.min()#最短距离 
    result_this.append(distance_all.min())#绘图用
    result_best.append(best)#绘图用

    #信息素更新
    tao=tao*0.5#挥发（数字越大，挥发越慢）
    ants_matrix1=np.concatenate((ants_matrix[:,1:],ants_matrix[:,:1]),axis=1)
    for i in range(city_num):
        tao[ants_matrix1[i],ants_matrix[i]]+=1#蚁数
    tao+=.001#人为给所有路段加上少量信息素，避免彻底收敛

print(r"最优解%.2f"%best,list(best_route))

#绘图
import matplotlib.pyplot as plt
plt.plot(range(Generations),result_best,alpha=0.5,color="green")#已知最优值曲线
plt.plot(range(Generations),result_this,alpha=0.5,color="red")#本代最优值曲线
plt.show()