有限级信息素蚁群算法使用路径等级作为信息素更新的依据，相比于传统的蚁群算法，舍弃了目标函数作为信息素更新的依据这一方法。在TSP问题中，目标函数实际就是路径长度，在传统的蚁群算法中信息素更新量为Q/f(x)，其中Q为某一常数，而f(x)就是目标函数值，即路径长度。而在有限级信息素蚁群算法中在这一点做出了改变。

首先给出算法的基本步骤：

步骤1：设定初始参数，初始化信息素

步骤2：按照路径选择规则构造问题的解

步骤3：按照信息素更新规则更新信息素

步骤4：判断是否达到停止条件，若满足，则停止迭代，否则跳至步骤2

有限级信息素蚁群算法主要变动就在步骤3，记为弧的级别，是一个单调递增函数，用于将弧的级别映射为信息素。在有限级信息素蚁群算法中，弧级别越大则弧上的信息素也就越高。表示路径惩罚等级，表示路径奖励等级。假设路径最小等级为1，路径最大等级为，为当前最优解，为当前次迭代最优解。

步骤3可以细分为以下几步：

<1>.

<2>.如果，则

<3>.对于，，即将当前最优解路径上的所有弧等级加，

<4>.，保证弧等级最小为1

<5>.，保证弧等级最大为

<6>.，通常选取为或者，为每条弧上的最大信息素

有限级信息素蚁群算法代码如下：

# coding: utf-8

# In[8]:


import numpy as np
from numpy import random as rd
import os
import matplotlib.pyplot as plt
np.set_printoptions(linewidth=1000, suppress=True)


# In[2]:

# 数据下载地址：https://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/tsp/
# 将数据存放在在代码存放目录的tspdata文件夹内
with open("./tspdata/kroA100.tsp", "r", encoding="utf-8") as file:
    data = file.read()
data = np.array([p.split(" ")[1:] for p in data.split("\n")[6:-2]]).astype(np.float32)
data


# In[3]:


def calc_distance_matrix(city_position_matrix):
    """
        city_position_matrix:城市坐标矩阵，形状为N行2列，N为城市数
        返回城市间距离矩阵distance_matrix，假设城市i到城市j的距离等于城市j到城市i的距离，因此，distance_matrix为上三角矩阵
    """
    city_number = city_position_matrix.shape[0]
    distance_matrix = np.zeros(shape=[city_number, city_number], dtype=np.float32)
    for i in range(city_number):
        from_city_position = city_position_matrix[i]
        for j in range(i + 1, city_number):
            to_city_position = city_position_matrix[j]
            distance_matrix[i, j] = np.sqrt(np.sum(np.power(city_position_matrix[i] - city_position_matrix[j], 2)))
    return distance_matrix


# In[4]:


class FGP_ACO(object):
    
    def __init__(self, initial_pheromone, r_2, r_1, M, beta, max_pheromone, run_times, city_position_matrix, ant_count, initial_route_level):
        """
            initial_pheromone:每条路径上的信息素初始值
            r_2:路径奖励等级
            r_1:路径惩罚级别
            M:路径最大级别
            beta:启发信息重要程度
            max_pheromone:信息素最大值
            run_times:最大循环次数
            city_position_matrix:城市坐标矩阵，形状为N行2列，N为城市数
            ant_count:蚂蚁数量
            initial_route_leve:表示每条路径的初始级别
        """
        self.r_2 = r_2
        self.r_1 = r_1
        self.M = M
        self.beta = beta
        self.max_pheromone = max_pheromone
        self.run_times = run_times
        self.city_position_matrix = city_position_matrix
        self.ant_count = ant_count
        # 城市间距离矩阵
        self.distance_matrix = calc_distance_matrix(city_position_matrix)
        # 城市数量
        self.N = city_position_matrix.shape[0]
        # 启发信息矩阵
        self.heuristic_information_matrix = np.triu(1 / self.distance_matrix, k=1)
        self.pheromone_matrix = np.triu(initial_pheromone * np.ones_like(self.distance_matrix), k=1)
        # 路径级别矩阵
        self.route_level_matrix = np.triu(np.ones_like(self.distance_matrix) * initial_route_level, k=1)
        # 每一代的最优距离列表
        self.optimal_distance_of_every_iteration = []
        # 每一代的当前最优距离列表
        self.current_optimal_distance_of_every_iteration = []
    
    def calc_move_prob_matrix(self, passed_city, not_passed_city):
        """
            passed_city:每只蚂蚁已经经过的城市
            not_passed_city:每只蚂蚁未经过的城市
            返回未经过城市中每个城市被选择的概率矩阵move_prob_matrix，形状为[self.ant_count, np.array(not_passed_city).shape[1]]
        """
        move_prob_matrix = []
        for i in range(self.ant_count):
            # 由于都是上三角矩阵，因此，如果路径为i到j，且i > j的，需要对掉，变为j到i
            start = np.array([passed_city[i][-1]] * len(not_passed_city[i]))
            end = np.array(not_passed_city[i])
            bool_index = start > end
            start[bool_index] = end[bool_index]
            end[bool_index] = passed_city[i][-1]
            # 找出每条未经过路径上的信息素not_passed_route_pheromone以及启发信息not_passed_route_heuristic_information，
            # 其长度均为未经过城市的个数
            not_passed_route_pheromone = self.pheromone_matrix[start, end]
            not_passed_route_heuristic_information = self.heuristic_information_matrix[start, end]
            move_prob_matrix.append(not_passed_route_pheromone * not_passed_route_heuristic_information ** self.beta / np.sum(not_passed_route_pheromone * not_passed_route_heuristic_information ** self.beta))
        return np.array(move_prob_matrix)
    
    def calc_route_distance(self, route):
        """
            route:表示一次循环的路径，例如[0, 2, 1, 4, 3]，表示0->2->1->4->3->0
            返回这条路径的长度
        """
        distance = 0
        for i in range(len(route) - 1):
            min_value, max_value = sorted([route[i], route[i + 1]])
            distance += self.distance_matrix[min_value, max_value]
        min_value, max_value = sorted([route[0], route[1]])
        distance += self.distance_matrix[min_value, max_value]
        return distance
    
    def g_function(self):
        """
            更新信息素矩阵
        """
        self.pheromone_matrix = np.sqrt((self.max_pheromone ** 2 - 1) / (self.M - 1) * (self.route_level_matrix - 1) + 1)
        
    def run(self):
        """
            迭代寻优，返回最优路径以及最优路径长度
        """
        # 随即初始化一个当前最优解, current_optimal_route表示当前最优路径，current_optimal_distance表示当前最优路径的长度
        current_optimal_route = rd.permutation(np.arange(self.N))
        current_optimal_distance = self.calc_route_distance(current_optimal_route)
        self.current_optimal_distance_of_every_iteration.append(current_optimal_distance)
        # 控制循环代数
        for times in range(1, 1 + self.run_times):
            if times % int(0.1 * self.run_times) == 0:
                print("第%d次迭代：" % times)
                print("最优路径：", current_optimal_route)
                print("最优路径长度：", current_optimal_distance)
            # 每只蚂蚁已经经过的城市列表[lst1, lst2, lst3, ...]，lst1表示第一只蚂蚁经过城市列表，lst1[-1]表示蚂蚁1当前所在城市，以此类推
            passed_city = [[rd.randint(0, self.N)] for i in range(self.ant_count)]
            # 每只蚂蚁未经过的城市列表[lst1, lst2, ...]，lst1表示第一只蚂蚁未经过的城市列表，以此类推
            not_passed_city = [list(set(range(self.N)) - set(i)) for i in passed_city]
            # 判断每只蚂蚁是否均遍历完所有城市，如果遍历完则跳出循环
            while len(np.unique(not_passed_city)) != 0:
                # 每次循环所有蚂蚁选择一个城市进行状态转移
                # 计算未经过城市中每个城市被选择的概率矩阵move_prob_matrix，形状为[self.ant_count, np.array(not_passed_city).shape[1]]
                move_prob_matrix = self.calc_move_prob_matrix(passed_city, not_passed_city)
                # 求取被选择概率最大的城市索引
                select_city_index = np.argmax(move_prob_matrix, axis=1)
                # 进行状态转移，更新已经经过的城市列表以及未经过的城市列表
                for i in range(self.ant_count):
                    passed_city[i].append(not_passed_city[i].pop(select_city_index[i]))
            # 更新路径等级，首先进行惩罚
            self.route_level_matrix -= self.r_1
            self.route_level_matrix[self.route_level_matrix < 1] = 1
            # 更新当前最优路径以及当前最优路径长度
            # 计算所有蚂蚁经过的路径的路径长度
            distance_lst = [self.calc_route_distance(r) for r in passed_city]
            # 最优蚂蚁索引
            optimal_index = np.argmin(distance_lst)
            current_iter_optimal_route = passed_city[optimal_index]
            current_iter_optimal_distance = distance_lst[optimal_index]
            self.optimal_distance_of_every_iteration.append(current_iter_optimal_distance)
            if current_iter_optimal_distance < current_optimal_distance:
                current_optimal_distance = current_iter_optimal_distance
                current_optimal_route = current_iter_optimal_route
            self.current_optimal_distance_of_every_iteration.append(current_optimal_distance)
            # 更新路径等级，进行奖励
            for i in range(len(current_optimal_route) - 1):
                min_value, max_value = sorted([current_optimal_route[i], current_optimal_route[i + 1]])
                self.route_level_matrix[min_value, max_value] += self.r_2
            min_value, max_value = sorted([current_optimal_route[0], current_optimal_route[-1]])
            self.route_level_matrix[min_value, max_value] += self.r_2
            self.route_level_matrix[self.route_level_matrix > self.M] = self.M
            # 更新信息素
            self.g_function()
#             print("最优路径：", current_optimal_route)
#             print("最优路径长度：", current_optimal_distance)
        return current_optimal_route, current_optimal_distance


# In[6]:


fgp_aco = FGP_ACO(100, 3, 1, 50, 5, 50, 2000, data, 25, 100)
print("城市坐标：\n", fgp_aco.city_position_matrix)
print("城市间距离矩阵：\n", fgp_aco.distance_matrix)
print("信息素矩阵：\n", fgp_aco.pheromone_matrix)
print("启发信息矩阵：\n", fgp_aco.heuristic_information_matrix)
fgp_aco.run()


# In[16]:


ax = plt.subplot(1, 1, 1)
ax.scatter(np.arange(len(fgp_aco.optimal_distance_of_every_iteration)), fgp_aco.optimal_distance_of_every_iteration, alpha=0.3)
plt.show()