有限級資訊素蟻群演算法使用路徑等級作為資訊素更新的依據，相比於傳統的蟻群演算法，捨棄了目標函式作為資訊素更新的依據這一方法。在TSP問題中，目標函式實際就是路徑長度，在傳統的蟻群演算法中資訊素更新量為Q/f(x)，其中Q為某一常數，而f(x)就是目標函式值，即路徑長度。而在有限級資訊素蟻群演算法中在這一點做出了改變。

首先給出演算法的基本步驟：

步驟1：設定初始引數，初始化資訊素

步驟2：按照路徑選擇規則構造問題的解

步驟3：按照資訊素更新規則更新資訊素

步驟4：判斷是否達到停止條件，若滿足，則停止迭代，否則跳至步驟2

有限級資訊素蟻群演算法主要變動就在步驟3，記為弧的級別，

步驟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()