Python实现蚁群算法

来源：脚本之家时间：2022-03-04 16:50:41

1、引言

在自然界中各种生物群体显现出来的智能近几十年来得到了学者们的广泛关注，学者们通过对简单生物体的群体行为进行模拟，进而提出了群智能算法。其中，模拟蚁群觅食过程的蚁群优化算法(Ant Colony Optimization, ACO)和模拟鸟群运动方式的粒子群算法(Particle Swarm Optimization， PSO)是两种最主要的群智能算法。

蚁群算法是一种源于大自然生物世界的新的仿生进化算法，由意大利学者M. Dorigo, V. Maniezzo和A.Colorni等人于20世纪90年代初期通过模拟自然界中蚂蚁集体寻径行为而提出的一种基于种群的启发式随机搜索算法"。蚂蚁有能力在没有任何提示的情形下找到从巢穴到食物源的最短路径，并且能随环境的变化，适应性地搜索新的路径，产生新的选择。其根本原因是蚂蚁在寻找食物时，能在其走过的路径上释放一种特殊的分泌物——信息素(也称外激素)，随着时间的推移该物质会逐渐挥发，后来的蚂蚁选择该路径的概率与当时这条路径上信息素的强度成正比。当一条路径上通过的蚂蚁越来越多时，其留下的信息素也越来越多，后来蚂蚁选择该路径的概率也就越高，从而更增加了该路径上的信息素强度。而强度大的信息素会吸引更多的蚂蚁，从而形成一种正反馈机制。通过这种正反馈机制，蚂蚁最终可以发现最短路径。

最早的蚁群算法是蚂蚁系统(Ant System,AS)，研究者们根据不同的改进策略对蚂蚁系统进行改进并开发了不同版本的蚁群算法，并成功地应用于优化领域。用该方法求解旅行商(TSP)问题、分配问题、车间作业调度(job-shop)问题，取得了较好的试验结果。蚁群算法具有分布式计算、无中心控制和分布式个体之间间接通信等特征，易于与其他优化算法相结合，它通过简单个体之间的协作表现出了求解复杂问题的能力，已被广泛应用于求解优化问题。蚁群算法相对而言易于实现，且算法中并不涉及复杂的数学操作，其处理过程对计算机的软硬件要求也不高，因此对它的研究在理论和实践中都具有重要的意义。

目前，国内外的许多研究者和研究机构都开展了对蚁群算法理论和应用的研究，蚁群算法已成为国际计算智能领域关注的热点课题。虽然目前蚁群算法没有形成严格的理论基础，但其作为一种新兴的进化算法已在智能优化等领域表现出了强大的生命力。

2 蚁群算法理论

蚁群算法是对自然界蚂蚁的寻径方式进行模拟而得出的一种仿生算法。蚂蚁在运动过程中，能够在它所经过的路径上留下信息素进行信息传递，而且蚂蚁在运动过程中能够感知这种物质，并以此来指导自己的运动方向。因此，由大量蚂蚁组成的蚁群的集体行为便表现出一种信息正反馈现象:某一路径上走过的蚂蚁越多，则后来者选择该路径的概率就越大。

3 算法理论图解

(1)在自然界中，蚂蚁的食物源总是随机分散于蚁巢周围，在蚁群协调、分工、合作后总能找到一条通往食物源的最短路径。现实中，我们能观察到大量蚂蚁在巢穴与食物源之间形成近乎直线的路径，而不是曲线、圆等其他形状，如图（a）。

(2)蚂蚁群体不仅能完成复杂的任务，并且还能适应环境的变化，如在蚁群运动路线上突然出现障碍物时，一开始各只蚂蚁分布是均匀的，不管路径长短，蚂蚁总是先按同等概率选择各条路径，如图（b）。

(3)蚂蚁在运动过程中，能在其经过的路径上留下信息素，并且能感知到这种物质的存在及其强度，并以此指导自己运动的方向，蚂蚁倾向于信息素浓度高的方向移动。在相同时间内较短路径上的信息素量就遗留得较多，则选择较短路径得蚂蚁也随即增多，如图（c）。

(4)不难看出，由于大量蚂蚁组成得蚁群集体行为表现出的一种信息正反馈现象，在某一路径上走过的蚂蚁越多，则后来者选择该路径的概率就越大，蚂蚁个体质检就是通过这种信息交流机制来搜索食物，并最终沿着最短路径行进，如图（d）。

4 人工蚁群优化过程

基于以上真实蚁群寻找食物时的最优路径选择问题，可以构造人工蚁群，来解决最优化问题，如TSP问题。人工蚁群中把具有简单功能的工作单元看作蚂蚁。二者的相似之处在于都是优先选择信息素浓度大的路径。较短路径的信息素浓度高，所以能够最终被所有蚂蚁选择，也就是最终的优化结果。两者的区别在于人工蚁群有一定的记忆能力，能够记忆已经访问过的节点。同时，人工蚁群再选择下一条路径的时候是按一定算法规律有意识地寻找最短路径，而不是盲目的。例如在TSP问题中，可以预先知道当前城市到下一个目的地的距离。

在TSP问题的人工蚁群算法中，假设m只蚂蚁在图的相邻节点间移动，从而协作异步地得到问题的解。每只蚂蚁的一步转移概率由图中的每条边上的两类参数决定:一是信息素值，也称信息素痕迹;二是可见度，即先验值。

信息素的更新方式有两种:一是挥发，也就是所有路径上的信息素以一定的比率减少，模拟自然蚁群的信息素随时间挥发的过程;二是增强，给评价值“好”(有蚂蚁走过)的边增加信息素。

蚂蚁向下一个目标的运动是通过一个随机原则来实现的，也就是运用当前所在节点存储的信息，计算出下一步可达节点的概率，并按此概率实现一步移动，如此往复，越来越接近最优解。

蚂蚁在寻找过程中，或在找到一个解后，会评估该解或解的一部分的优化程度，并把评价信息保存在相关连接的信息素中。

这种算法有别于传统编程模式，其优势在于，避免了冗长的编程和筹划，程序本身是基于一定规则的随机运行来寻找最佳配置。也就是说，当程序最开始找到目标的时候，路径可能不是最优的。但是，程序可以通过蚂蚁寻找食物的时候的信息素原理，不断地去修正原来的路线，使整个路线越来越短，也就是说，程序执行的时间越长（在程序中也就是迭代次数不能太少，同时还要保证一定的蚂蚁数量），所获得的路径就越可能接近最优路径。这看起来很类似与我们所见的由无数例子进行归纳概括形成最佳路径的过程。实际上好似是程序的一个自我学习的过程。

这种优化过程的本质在于：

选择机制：信息素越多的路径，被选择的概率越大。

更新机制：路径上面的信息素会随蚂蚁的经过而增长，而且同时也随时间的推移逐渐挥发消失。

协调机制：蚂蚁间实际上是通过分泌物来互相通信、协同工作的。

蚁群算法正是充分利用了选择、更新和协调的优化机制，即通过个体之间的信息交流与相互协作最终找到最优解，使它具有很强的发现较优解的能力。

事实上，每只蚂蚁并不是像我们想象的需要知道整个世界的信息，他们其实只关心很小范围内的眼前信息，而且根据这些局部信息利用几条简单的规则进行决策，但是，当集群里有无数蚂蚁的时候，复杂性的行为就会凸现出来。这就是人工生命、复杂性科学解释的规律！那么，这些简单规则是什么呢？下面详细说明：

1、范围：

蚂蚁观察到的范围是一个方格世界，蚂蚁有一个参数为速度半径（一般是3），那么它能观察到的范围就是3*3个方格世界，并且能移动的距离也在这个范围之内。

2、环境：

蚂蚁所在的环境是一个虚拟的世界，其中有障碍物，有别的蚂蚁，还有信息素，信息素有两种，一种是找到食物的蚂蚁洒下的食物信息素，一种是找到窝的蚂蚁洒下的窝的信息素。每个蚂蚁都仅仅能感知它范围内的环境信息。环境以一定的速率让信息素消失。

3、觅食规则：

在每只蚂蚁能感知的范围内寻找是否有食物，如果有就直接过去。否则看是否有信息素，并且比较在能感知的范围内哪一点的信息素最多，这样，它就朝信息素多的地方走，并且每只蚂蚁多会以小概率犯错误，从而并不是往信息素最多的点移动。蚂蚁找窝的规则和上面一样，只不过它对窝的信息素做出反应，而对食物信息素没反应。

4、移动规则：

每只蚂蚁都朝向信息素最多的方向移，并且，当周围没有信息素指引的时候，蚂蚁会按照自己原来运动的方向惯性的运动下去，并且，在运动的方向有一个随机的小的扰动。为了防止蚂蚁原地转圈，它会记住最近刚走过了哪些点，如果发现要走的下一点已经在最近走过了，它就会尽量避开。

5、避障规则：

如果蚂蚁要移动的方向有障碍物挡住，它会随机的选择另一个方向，并且有信息素指引的话，它会按照觅食的规则行为。

6、播撒信息素规则：

每只蚂蚁在刚找到食物或者窝的时候撒发的信息素最多，并随着它走远的距离，播撒的信息素越来越少。

根据这几条规则，蚂蚁之间并没有直接的关系，但是每只蚂蚁都和环境发生交互，而通过信息素这个纽带，实际上把各个蚂蚁之间关联起来了。比如，当一只蚂蚁找到了食物，它并没有直接告诉其它蚂蚁这儿有食物，而是向环境播撒信息素，当其它的蚂蚁经过它附近的时候，就会感觉到信息素的存在，进而根据信息素的指引找到了食物。

那么，蚂蚁究竟是怎么找到食物的呢？

在没有蚂蚁找到食物的时候，环境没有有用的信息素，那么蚂蚁为什么会相对有效的找到食物呢？这要归功于蚂蚁的移动规则，尤其是在没有信息素时候的移动规则。首先，它要能尽量保持某种惯性，这样使得蚂蚁尽量向前方移动（开始，这个前方是随机固定的一个方向），而不是原地无谓的打转或者震动；其次，蚂蚁要有一定的随机性，虽然有了固定的方向，但它也不能像粒子一样直线运动下去，而是有一个随机的干扰。这样就使得蚂蚁运动起来具有了一定的目的性，尽量保持原来的方向，但又有新的试探，尤其当碰到障碍物的时候它会立即改变方向，这可以看成一种选择的过程，也就是环境的障碍物让蚂蚁的某个方向正确，而其他方向则不对。这就解释了为什么单个蚂蚁在复杂的诸如迷宫的地图中仍然能找到隐蔽得很好的食物。当然，在有一只蚂蚁找到了食物的时候，其他蚂蚁会沿着信息素很快找到食物的。

蚂蚁如何找到最短路径的？

这一是要归功于信息素，另外要归功于环境，具体说是计算机时钟。信息素多的地方显然经过这里的蚂蚁会多，因而会有更多的蚂蚁聚集过来。假设有两条路从窝通向食物，开始的时候，走这两条路的蚂蚁数量同样多（或者较长的路上蚂蚁多，这也无关紧要）。当蚂蚁沿着一条路到达终点以后会马上返回来，这样，短的路蚂蚁来回一次的时间就短，这也意味着重复的频率就快，因而在单位时间里走过的蚂蚁数目就多，洒下的信息素自然也会多，自然会有更多的蚂蚁被吸引过来，从而洒下更多的信息素……；而长的路正相反，因此，越来越多地蚂蚁聚集到较短的路径上来，最短的路径就近似找到了。也许有人会问局部最短路径和全局最短路的问题，实际上蚂蚁逐渐接近全局最短路的，为什么呢？这源于蚂蚁会犯错误，也就是它会按照一定的概率不往信息素高的地方走而另辟蹊径，这可以理解为一种创新，这种创新如果能缩短路途，那么根据刚才叙述的原理，更多的蚂蚁会被吸引过来。

5 基本蚁群算法及其流程

5.1蚁群算法公式

蚁群算法实际上是正反馈原理和启发式算法相结合的一种算法。在选择路径时，蚂蚁不仅利用了路径上的信息素，而且用到了城市间距离的倒数作为启发式因子。实验结果表明，ant-cycle模型比ant-quantity和ant-density模型有更好的性能。这是因为ant-cycle模型利用全局信息更新路径上的信息素量,而ant-quantity和ant-density模型使用局部信息。

5.2 蚁群算法程序概括

(1)参数初始化

在寻最短路钱，需对程序各个参数进行初始化，蚁群规模m、信息素重要程度因子α、启发函数重要程度因子β、信息素会发因子、最大迭代次数ddcs_max，初始迭代值为ddcs=1。

(2)构建解空间

将每只蚂蚁随机放置在不同的出发地点，对蚂蚁访问行为按照公式计算下一个访问的地点，直到所有蚂蚁访问完所有地点。

(3)更新信息素

计算每只蚂蚁经过的路径总长Lk,记录当前循环中的最优路径，同时根据公式对各个地点间连接路径上的信息素浓度进行更新。

(4)判断终止

迭代次数达到最大值前，清空蚂蚁经过的记录，并返回步骤2。

5.3 流程图

6 案例实现

6.1 案例1

求解函数：的最小值，其中x的取值范围为[-5，5], y的取值范围为[-5, 5]。这是一个有多个局部极值的函数。

6.2 Python实现

import numpy as np
from tqdm import tqdm#进度条设置
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib; matplotlib.use("TkAgg")
mpl.rcParams["font.sans-serif"] = ["SimHei"]  # 指定默认字体
mpl.rcParams["axes.unicode_minus"] = False  # 解决保存图像是负号"-"显示为方块的问题
 
 
#============蚁群算法求函数极值================
 
#=======适应度函数=====
def func(x,y):
    value = 20*np.power(x*x-y*y,2)-np.power(1-y,2)-3*np.power(1+y,2)+0.3
    return value
#=======初始化参数====
m=20                   #蚂蚁个数
G_max=200              #最大迭代次数
Rho=0.9                #信息素蒸发系数
P0=0.2                 #转移概率常数
XMAX= 5               #搜索变量x最大值
XMIN= -5              #搜索变量x最小值
YMAX= 5               #搜索变量y最大值
YMIN= -5              #搜索变量y最小值
X=np.zeros(shape=(m,2)) #蚁群 shape=(20, 2)
Tau=np.zeros(shape=(m,)) #信息素
P=np.zeros(shape=(G_max,m)) #状态转移矩阵
fitneess_value_list=[] #迭代记录最优目标函数值
#==随机设置蚂蚁初始位置==
for i in range(m):#遍历每一个蚂蚁
    X[i,0]=np.random.uniform(XMIN,XMAX,1)[0] #初始化x
    X[i,1]=np.random.uniform(YMIN,YMAX,1)[0] #初始化y
    Tau[i]=func(X[i,0],X[i,1])
 
step=0.1;                #局部搜索步长
for NC in range(G_max):#遍历每一代
    lamda=1/(NC+1)
    BestIndex=np.argmin(Tau) #最优索引
    Tau_best=Tau[BestIndex] #最优信息素
     #==计算状态转移概率===
    for i in range(m):#遍历每一个蚂蚁
        P[NC,i]=np.abs((Tau_best-Tau[i]))/np.abs(Tau_best)+0.01 #即例最优信息素的距离
 
    #=======位置更新==========
    for i in range(m):  # 遍历每一个蚂蚁
        #===局部搜索====
        if P[NC,i] XMAX:
            temp1 =XMAX
        if temp2 < XMIN:
            temp2 =XMIN
        if temp2 > XMAX:
            temp2 =XMAX
 
 
        #==判断蚂蚁是否移动(选更优===
        if func(temp1, temp2) < func(X[i, 0], X[i, 1]):
            X[i, 0] = temp1
            X[i, 1]= temp2
 
    #=====更新信息素========
    for i in range(m):  # 遍历每一个蚂蚁
        Tau[i] = (1 - Rho) * Tau[i] + func(X[i, 0], X[i, 1]) #(1 - Rho) * Tau[i] 信息蒸发后保留的
        index=np.argmin(Tau)#最小值索引
        value=Tau[index]#最小值
    fitneess_value_list.append(func(X[index,0],X[index,1])) #记录最优目标函数值
 
 
 
#==打印结果===
min_index=np.argmin(Tau)#最优值索引
minX=X[min_index,0]  #最优变量x
minY=X[min_index,1]  #最优变量y
minValue=func(X[min_index,0],X[min_index,1])  #最优目标函数值
 
print("最优变量x",minX,end="")
print("最优变量y",minY,end="\n")
print("最优目标函数值",minValue)
 
plt.plot(fitneess_value_list,label="迭代曲线")
plt.legend()
plt.show()

6.3 结果

最优变量x 5.0最优变量y 5.0
最优目标函数值 -123.7

6.4 案例2

6.5 Python实现

#====================导入相关库=============================
import pandas as pd
import numpy as np
from tqdm import tqdm#进度条设置
import matplotlib.pyplot as plt
import matplotlib; matplotlib.use("TkAgg")
from pylab import *
mpl.rcParams["font.sans-serif"] = ["SimHei"]
mpl.rcParams["axes.unicode_minus"] = False
 
 
#=======================定义函数==========================
 
#=======目标函数=====
def calc_f(X):
    """计算粒子的的适应度值，也就是目标函数值，X 的维度是 size * 2 """
    A = 10
    pi = np.pi
    x = X[0]
    y = X[1]
    return 2 * A + x ** 2 - A * np.cos(2 * pi * x) + y ** 2 - A * np.cos(2 * pi * y)
 
#====惩罚项函数======
def calc_e(X):
    """计算蚂蚁的惩罚项，X 的维度是 size * 2 """
    ee = 0
    """计算第一个约束的惩罚项"""
    e1 = X[0] + X[1] - 6
    ee += max(0, e1)
    """计算第二个约束的惩罚项"""
    e2 = 3 * X[0] - 2 * X[1] - 5
    ee += max(0, e2)
    return ee
 
#===子代和父辈之间的选择操作====
def update_best(parent,parent_fitness,parent_e,child,child_fitness,child_e):
    """
        针对不同问题，合理选择惩罚项的阈值。本例中阈值为0.00001
        :param parent: 父辈个体
        :param parent_fitness:父辈适应度值
        :param parent_e    ：父辈惩罚项
        :param child:  子代个体
        :param child_fitness 子代适应度值
        :param child_e  ：子代惩罚项
        :return: 父辈 和子代中较优者、适应度、惩罚项
        """
    # 规则1，如果 parent 和 child 都没有违反约束，则取适应度小的
    if parent_e <= 0.00001 and child_e <= 0.00001:
        if parent_fitness <= child_fitness:
            return parent,parent_fitness,parent_e
        else:
            return child,child_fitness,child_e
    # 规则2，如果child违反约束而parent没有违反约束，则取parent
    if parent_e < 0.00001 and child_e  >= 0.00001:
        return parent,parent_fitness,parent_e
    # 规则3，如果parent违反约束而child没有违反约束，则取child
    if parent_e >= 0.00001 and child_e < 0.00001:
        return child,child_fitness,child_e
    # 规则4，如果两个都违反约束，则取适应度值小的
    if parent_fitness <= child_fitness:
        return parent,parent_fitness,parent_e
    else:
        return child,child_fitness,child_e
 
#=======================初始化参数==========================
m=20                   #蚂蚁个数
G_max=200              #最大迭代次数
Rho=0.9                #信息素蒸发系数
P0=0.2                 #转移概率常数
XMAX= 2               #搜索变量x最大值
XMIN= 1              #搜索变量x最小值
YMAX= 0               #搜索变量y最大值
YMIN= -1              #搜索变量y最小值
step=0.1                #局部搜索步长
P=np.zeros(shape=(G_max,m)) #状态转移矩阵
fitneess_value_list=[] #迭代记录最优目标函数值
 
#=======================初始化蚂蚁群体位置和信息素==========================
def initialization():
    """
    :return: 初始化蚁群和初始信息素
    """
    X = np.zeros(shape=(m, 2))  # 蚁群 shape=(20, 2)
    Tau = np.zeros(shape=(m,))  # 信息素
    for i in range(m):  # 遍历每一个蚂蚁
        X[i, 0] = np.random.uniform(XMIN, XMAX, 1)[0]  # 初始化x
        X[i, 1] =np.random.uniform(YMIN, YMAX, 1)[0]  # 初始化y
        Tau[i] = calc_f(X[i])#计算信息素
    return X,Tau
 
#===位置更新====
def position_update(NC,P,X):
    """
    :param NC: 当前迭代次数
    :param P: 状态转移矩阵
    :param X: 蚁群
    :return: 蚁群X
    """
    lamda = 1 / (NC + 1)
    # =======位置更新==========
    for i in range(m):  # 遍历每一个蚂蚁
        # ===局部搜索===
        if P[NC, i] < P0:
            temp1 = X[i, 0] + (2 * np.random.random() - 1) * step * lamda  # x(2 * np.random.random() - 1) 转换到【-1,1】区间
            temp2 = X[i, 1] + (2 * np.random.random() - 1) * step * lamda  # y
        # ===全局搜索===
        else:
            temp1 = X[i, 0] + (XMAX - XMIN) * (np.random.random() - 0.5)
            temp2 = X[i, 0] + (YMAX - YMIN) * (np.random.random() - 0.5)
 
        # =====边界处理=====
        if (temp1 < XMIN) or (temp1 > XMAX):
            temp1 = np.random.uniform(XMIN, XMAX, 1)[0]  # 初始化x
        if (temp2 < YMIN) or (temp2 > YMAX):
            temp2 = np.random.uniform(YMIN, YMAX, 1)[0]  # 初始化y
 
        #=====判断蚂蚁是否移动(选更优）=====
        #==子代蚂蚁==
        children=np.array([temp1,temp2])#子代个体蚂蚁
        children_fit=calc_f(children) #子代目标函数值
        children_e=calc_e(children) #子代惩罚项
        parent=X[i]#父辈个体蚂蚁
        parent_fit=calc_f(parent)#父辈目标函数值
        parent_e=calc_e(parent)#父辈惩罚项
 
        pbesti, pbest_fitness, pbest_e = update_best(parent, parent_fit, parent_e, children, children_fit,children_e)
        X[i]=pbesti
    return X
 
#======信息素更新============
def Update_information(Tau,X):
    """
    :param Tau: 信息素
    :param X: 蚂蚁群
    :return: Tau信息素
    """
    for i in range(m):  # 遍历每一个蚂蚁
        Tau[i] = (1 - Rho) * Tau[i] + calc_f(X[i]) #(1 - Rho) * Tau[i] 信息蒸发后保留的
    return Tau
 
#=============主函数======================
def main():
    X,Tau=initialization() #初始化蚂蚁群X 和信息素 Tau
    for NC in tqdm(range(G_max)):  # 遍历每一代
        BestIndex = np.argmin(Tau)  # 最优索引
        Tau_best = Tau[BestIndex]  # 最优信息素
        # 计算状态转移概率
        for i in range(m):  # 遍历每一个蚂蚁
            P[NC, i] = np.abs((Tau_best - Tau[i])) / np.abs(Tau_best) + 0.01  # 即离最优信息素的距离
        # =======位置更新==========
        X=position_update(NC,P,X) #X.shape=(20, 2)
 
        # =====更新信息素========
        Tau=Update_information(Tau, X)
 
        # =====记录最优目标函数值========
        index = np.argmin(Tau)  # 最小值索引
        value = Tau[index]  # 最小值
        fitneess_value_list.append(calc_f(X[index]))  # 记录最优目标函数值
 
    #=====打印结果=======
    min_index = np.argmin(Tau)  # 最优值索引
    minX = X[min_index, 0]  # 最优变量x
    minY = X[min_index, 1]  # 最优变量y
    minValue = calc_f(X[min_index])  # 最优目标函数值
 
    print("最优变量x", minX, end="") 
    print("最优变量y", minY, end="\n")
    print("最优目标函数值", minValue)
    print("最优变量对应的惩罚项",calc_e(X[min_index]))
    #=====可视化=======
    plt.plot(fitneess_value_list, label="迭代曲线")
    plt.legend()
    plt.show()
 
 
 
if __name__=="__main__":
    main()