【数值优化】启发搜索之群体智能（2）-- 蚁群

蚁群算法（ACO）是一种模拟蚂蚁觅食行为的群体智能优化算法，由意大利学者Dorigo、Maniezzo等人于1991年提出，并首先使用在解决TSP（旅行商问题）上。
蚁群觅食过程中，每只蚂蚁在所走过的路径上均会释放出一种信息素，该信息素随时间的推移逐渐挥发。因此，每条路径上的信息素同时存在正负反馈两种机制。正反馈：蚂蚁每次经过该路径均会释放信息素使得该路径上的信息素浓度增加；负反馈：每条路径上的信息素随时间推移会逐渐挥发。由此，我们可以判断，在起点与终点之间，当相同数量的蚂蚁初始同时经过两条不同的路径时，路径上初始信息素的浓度是相同的；不过，当路径越短时，信息素挥发时间也越短，残留信息素浓度也将越高。随后的蚂蚁将根据路径上残留信息素浓度的大小对路径进行选择 — 浓度越高，选择概率越大。最终导致信息素浓度越高的路径上蚂蚁的选择数目越多，而更多的蚂蚁也将同时导致该路径上残留信息素浓度越高（即高者越高，低者越低）。因此，在理想情况下，整个蚁群将逐渐向信息素浓度最高的路径（即最短路径）进行转移。

蚁群算法过程

旅行商问题（travelling salesman problem，TSP）是物流领域的典型问题。蚂蚁算法求解TSP问题的过程如下：

STEP1： 随机产生种群，设迭代的最大次数为itermax，初始iter=0。
STEP2： 将N_ant只蚂蚁置于N_city个顶点上。
STEP3： N_ant只蚂蚁按概率函数选择下一个城市，完成各自的周游。
STEP4： 记录本次迭代最佳路线。
STEP5： 全局更新信息素值。
STEP6： 若终止条件满足，则结束；否则iter=iter+1，转入步骤2进行下一代优化。终止条件可以是迭代的代数，也可限定运行的时间，或设定最短路长的下限。
STEP7： 输出结果

蚁群算法的代码实现

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

#========================================================
#  数据准备
#========================================================

# 初始化城市坐标，总共52个城市
coordinates = np.array([[565.0, 575.0], [25.0, 185.0], [345.0, 750.0], [945.0, 685.0], [845.0, 655.0],
                        [880.0, 660.0], [25.0, 230.0], [525.0, 1000.0], [580.0, 1175.0], [650.0, 1130.0],
                        [1605.0, 620.0], [1220.0, 580.0], [1465.0, 200.0], [1530.0, 5.0], [845.0, 680.0],
                        [725.0, 370.0], [145.0, 665.0], [415.0, 635.0], [510.0, 875.0], [560.0, 365.0],
                        [300.0, 465.0], [520.0, 585.0], [480.0, 415.0], [835.0, 625.0], [975.0, 580.0],
                        [1215.0, 245.0], [1320.0, 315.0], [1250.0, 400.0], [660.0, 180.0], [410.0, 250.0],
                        [420.0, 555.0], [575.0, 665.0], [1150.0, 1160.0], [700.0, 580.0], [685.0, 595.0],
                        [685.0, 610.0], [770.0, 610.0], [795.0, 645.0], [720.0, 635.0], [760.0, 650.0],
                        [475.0, 960.0], [95.0, 260.0], [875.0, 920.0], [700.0, 500.0], [555.0, 815.0],
                        [830.0, 485.0], [1170.0, 65.0], [830.0, 610.0], [605.0, 625.0], [595.0, 360.0],
                        [1340.0, 725.0], [1740.0, 245.0]])

#========================================================
#  辅助工具
#========================================================

def getdistmat(coordinates):
    '''
    计算52个城市间的欧式距离矩阵
    '''
    num = coordinates.shape[0]
    distmat = np.zeros((52, 52)) # 初始化生成52*52的矩阵
    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

#========================================================
#  蚁群算法超参数
#========================================================
 
numant = 60  # 蚂蚁个数
numcity = coordinates.shape[0]  # 城市个数, 也就是任务个数
alpha = 1  # 信息素重要程度因子
beta = 5   # 启发函数重要程度因子
rho = 0.1  # 信息素的挥发速度
Q = 1      # 完成率

#========================================================
#  核心算法
#========================================================

class ACO:
    '''
    蚁群算法
    '''
    def __init__(self, distmat, itermax):
        '''
        distmat: 城市的距离矩阵
        itermax: 最大迭代次数
        '''
        self.itermax = itermax
        self.etatable = 1.0 / (distmat + np.diag([1e10] * numcity))  # 启发函数矩阵，表示蚂蚁从城市i转移到城市j的期望程度
        self.pheromonetable = np.ones((numcity, numcity))  # 信息素矩阵
        self.pathtable = np.zeros((numant, numcity)).astype(int) # 路径记录表，转化成整型
        self.lengthaver = np.zeros(self.itermax)  # 存放每次迭代后，路径的平均长度
        self.lengthbest = np.zeros(self.itermax)  # 存放每次迭代后，最佳路径长度
        self.pathbest = np.zeros((self.itermax, numcity))  # 存放每次迭代后，最佳路径城市的坐标

    def showhistory(self):
        '''
        绘制蚁群算法迭代过程
        '''
        fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(12, 10))
        axes[0].plot(self.lengthaver, 'k', marker='*')  # 平均路径长度
        axes[0].set_title('Average Length')  
        axes[0].set_xlabel(u'self.iteration')
        axes[1].plot(self.lengthbest, 'k', marker='<')  # 最优路径长度
        axes[1].set_title('Best Length')
        axes[1].set_xlabel(u'self.iteration')
        plt.show()
    
    def savebestpath(self):
        '''
        保存最优路径
        '''
        # 作出找到的最优路径图
        bestpath = self.pathbest[-1]
        plt.plot(coordinates[:, 0], coordinates[:, 1], 'r.', marker='>')
        plt.xlim([-100, 2000])
        plt.ylim([-100, 1500])
        for i in range(numcity - 1):
            # 按坐标绘出最佳两两城市间路径
            m, n = int(bestpath[i]), 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(bestpath[51])][0]],    [coordinates[int(bestpath[0])][1],coordinates[int(bestpath[50])][1]]   ,'b')
        
        ax = plt.gca()
        ax.set_title("Best Path")
        ax.set_xlabel('X_axis')
        ax.set_ylabel('Y_axis')
        
        plt.savefig('Best Path.png', dpi=500, bbox_inches='tight')
        plt.close()

    def main(self):
        '''
        主流程:蚁群迭代过程
        '''
        self.iter = 0
        while self.iter < self.itermax:
            # 将蚂蚁随机放置于各个城市中
            if numant <= numcity:  # 城市数比蚂蚁数多，不用管
                self.pathtable[:, 0] = np.random.permutation(range(numcity))[:numant]  # 矩阵表示哪个蚂蚁在哪个城市
            else:  # 蚂蚁数比城市数多，需要有城市放多个蚂蚁
                self.pathtable[:numcity, 0] = np.random.permutation(range(numcity))[:]
                self.pathtable[numcity:, 0] = np.random.permutation(range(numcity))[:numant - numcity]
            
            # 算出第i只蚂蚁转移到下一个城市的概率
            length = np.zeros(numant)
            for i in range(numant):
                visiting = self.pathtable[i, 0]  # 当前所在的城市
                unvisited = set(range(numcity))  # 未访问的城市集合
                unvisited.remove(visiting)  # 删除已经访问过的城市元素
                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(self.pheromonetable[visiting][listunvisited[k]], alpha) * np.power(self.etatable[visiting][listunvisited[k]], alpha)
                    # 概率公式的分母,其中[visiting][listunvis[k]]是从本城市到k城市的信息素
                    cumsumprobtrans = (probtrans / sum(probtrans)).cumsum()
                    # 本只蚂蚁的转移到各个城市的概率斐波衲挈数列
                    cumsumprobtrans -= np.random.rand()
                    # 随机生成下个城市的转移概率，再用区间比较
                    k = listunvisited[list(cumsumprobtrans > 0).index(True)]
                    self.pathtable[i, j] = k  # 下一个要访问的城市
                    # 采用禁忌表来记录蚂蚁i当前走过的第j城市的坐标, 这里走了第j个城市. k是中间值
                    unvisited.remove(k)
                    # 将未访问城市列表中的K城市删去，增加到已访问城市列表中
                    length[i] += distmat[visiting][k]
                    # 计算本城市到K城市的距离
                    visiting = k
                # 计算本只蚂蚁的总的路径距离，包括最后一个城市和第一个城市的距离
                length[i] += distmat[visiting][self.pathtable[i, 0]]
            # 本轮的平均路径
            self.lengthaver[self.iter] = length.mean()
            
            # 求出最佳路径
            if self.iter == 0:  # 如果是第一轮路径，则选择本轮最短的路径并返回索引值下标, 将其记录
                self.lengthbest[self.iter] = length.min()
                self.pathbest[self.iter] = self.pathtable[length.argmin()].copy()
            else: # 后面几轮的情况，更新最佳路径
                if length.min() > self.lengthbest[self.iter - 1]:
                    self.lengthbest[self.iter] = self.lengthbest[self.iter - 1]
                    self.pathbest[self.iter] = self.pathbest[self.iter - 1].copy()
                else:
                    self.lengthbest[self.iter] = length.min()
                    self.pathbest[self.iter] = self.pathtable[length.argmin()].copy()
        
            # 更新信息素
            changepheromonetable = np.zeros((numcity, numcity))
            # 更新所有的蚂蚁
            for i in range(numant): 
                for j in range(numcity - 1):
                    # 根据公式更新本只蚂蚁改变的城市间的信息素
                    changepheromonetable[self.pathtable[i, j]][self.pathtable[i, j + 1]] += Q / distmat[self.pathtable[i, j]][self.pathtable[i, j + 1]]
                # 首城市到最后一个城市 所有蚂蚁改变的信息素总和
                changepheromonetable[self.pathtable[i, j + 1]][self.pathtable[i, 0]] += Q / distmat[self.pathtable[i, j + 1]][self.pathtable[i, 0]]
            # 信息素更新公式p=(1-挥发速率)*现有信息素+改变的信息素
            self.pheromonetable = (1 - rho) * self.pheromonetable + changepheromonetable
            self.iter += 1  # 迭代次数指示器+1
        # 打印迭代过程
        self.showhistory()
        # 保存最优路径
        self.savebestpath()

if __name__ == "__main__":
    distmat = getdistmat(coordinates) # 城市的距离矩阵
    aco = ACO(distmat, 100)
    aco.main()