任务要求

写论文和报告（论文流程：问题——算法——评估——结论)

实践

代码

solutionbase.py

#此为蚁群算法基线模型
import numpy as np
import random as ra

def TSP(net):
    '''初始化'''
    n = net.shape[0]
    h = 0.1 #挥发率
    m = 3 #蚂蚁数量
    t0 = m / TanXin(n, net) #初始信息素浓度
    a = 1
    b = 2
    #存储结构
    R = np.ones([m, n+1], dtype=int) #蚂蚁路径矩阵
    T= np.ones([n, n]) * t0   #信息浓度矩阵
    C = np.ones(m)        #路径长度序列
    global_best = float('inf') #全局最优值
    global_bestR = None #全局最优路径

    '''循环'''
    for turn in range(50):
        for k in range(m):  # 蚂蚁k
            for p in range(n + 1):
                if p == 0:
                    R[k][0] = ra.randint(0, n-1)
                elif p == n:
                    R[k][n] = R[k][0]
                else:
                    i = R[k][p - 1]  # 前一个结点
                    js = neighbor(n, net, i)
                    js = list(set(js) - set(R[k][:p]))  # 下一个可能节点
                    tjs = [T[i, j]**a * (1 / net[i, j])**b for j in js]
                    sum_tjs = sum(tjs)
                    pjs = [tj / sum_tjs for tj in tjs]  # 下一个可能节点的概率
                    # 轮盘算法求下一个节点j
                    r = ra.random()
                    for z in range(len(pjs)):
                        if r <= pjs[z]:
                            j = js[z]
                        else:
                            r -= pjs[z]
                    R[k][p] = j  # 把j加入k蚂蚁的路径
            C[k] = sum([net[R[k][i], R[k][i+1]] for i in range(n)])
        # 更新t
        for i in range(n):
            for j in range(n):
                T[i, j] = (1-h)* T[i, j]
                for k in range(m):
                    if j == R[k][list(R[k]).index(i)+1]:
                        T[i, j] += 1 / C[k]
        #更新全局最优值
        local_best = min(C) #局部最优值
        bestk = list(C).index(local_best) #局部最优值对应的蚂蚁编号
        if local_best < global_best:
            global_best = local_best
            global_bestR = R[bestk].copy()
            

    # '''输出'''
    # print('最短路径为：', global_best)
    # print('最短路径长度：', global_bestR)
    return global_best


def TanXin(n, net):
    R = np.ones(n+1, dtype=int)
    while True:
        for p in range(n+1):
            if p == 0:
                R[0] = ra.randint(0, n-1)
            elif p == n:
                R[n] = R[0]
            else:
                i = R[p-1]
                js = neighbor(n, net, i)
                js = list(set(js)-set(R[:p]))
                ljs = [net[i, j] for j in js]
                j = js[ljs.index(min(ljs))]
                R[p] = j
        C = sum([net[R[i], R[i+1]] for i in range(n)])
        if C != float('inf'):
            break
    return C

def neighbor(n, net, i):
    return [j for j in range(n) if j != i if net[i][j] != float('inf')]

if __name__ == '__main__':
    net = np.array([[float('inf'), 3, 1, 2],
                    [3, float('inf'), 5, 4],
                    [1, 5, float('inf'), 2],
                    [2, 4, 2, float('inf')]])
    TSP(net)

solution.py

#此为蚁群算法的改进模型
'''
改进:
(1) 信息素初始化按根据距离分配，加快收敛(具体使得平均值为t0,但是具体每一个值按路径长度分配)
(2) 设置稳定时刻点turn0，如果在超过规定时间(y)全局最优值没有发生变化，那么此时为turn0时刻
turn0时刻后的信息素挥发速率变快，蚁群产生信息素减少，以此增加路径选择的发散性
具体为：信息素保留率：(1-h) * 1 / (turn-turn0)
       一只蚂蚁在(i, j)上生产的信息素： 1 / C[k] * 1 / (turn-turn0)
(3) 设置重来时刻点turn1,如果在turn0后的一段时间(x)后全局最优值没有发生变化，那么此时为turn1时刻
此时重新初始化，融断最多发生u次
'''

import numpy as np
import random as ra

def TSP(net):
    '''初始化'''
    n = net.shape[0]
    h = 0.1 #挥发率
    m = 3 #蚂蚁数量
    t0 = m / TanXin(n, net) #初始信息素浓度
    a = 1
    b = 2
    #存储结构
    R = np.ones([m, n+1], dtype=int) #蚂蚁路径矩阵
    T = np.ones([n, n])  #信息浓度矩阵 
    #初始化信息浓度矩阵
    sum_edge = 0
    num_edge = 0
    for i in range(n):
        for j in range(n):
            if i == j or net[i, j] == float('inf'):
                continue
            sum_edge += 1.0 / net[i, j]
            num_edge += 1;
    for i in range(n):
        for j in range(n):
            if i == j or net[i, j] == float('inf'):
                continue
            T[i, j] = T[i, j] * num_edge * t0 * 1 / net[i, j] / sum_edge
    C = np.ones(m)        #路径长度序列
    global_best = float('inf') #全局最优值
    global_bestR = None #全局最优路径

    #改进部分 
    y = 5
    x = 5
    u = 5
    turn0times = 0  #计时器1（记录到turn0前未变化的次数）
    turn1times = 0  #计时器2 (记录到turn0后为变化的次数)
    turn1num = 0 #计数器2 (记录重来的次数)
    turn0 = None #稳定时刻点
    isturn0 = False #是否到达稳定时刻点
    isturn1 = True #是否到达重来时刻点
  


    '''循环'''
    turn = 0
    while turn <= 50:
        if isturn0:
            turn1times += 1
        for k in range(m):  # 蚂蚁k
            for p in range(n + 1):
                if p == 0:
                    R[k][0] = ra.randint(0, n-1)
                elif p == n:
                    R[k][n] = R[k][0]
                else:
                    i = R[k][p - 1]  # 前一个结点
                    js = neighbor(n, net, i)
                    js = list(set(js) - set(R[k][:p]))  # 下一个可能节点
                    tjs = [T[i, j]**a * (1 / net[i, j])**b for j in js]
                    sum_tjs = sum(tjs)
                    pjs = [tj / sum_tjs for tj in tjs]  # 下一个可能节点的概率
                    # 轮盘算法求下一个节点j
                    r = ra.random()
                    for z in range(len(pjs)):
                        if r <= pjs[z]:
                            j = js[z]
                        else:
                            r -= pjs[z]
                    R[k][p] = j  # 把j加入k蚂蚁的路径
            C[k] = sum([net[R[k][i], R[k][i+1]] for i in range(n)])
        #更新T
        for i in range(n):
            for j in range(n):
                T[i, j] = (1-h)* T[i, j]
                if isturn0:
                    T[i, j] *= 1.0 / (turn-turn0)
                for k in range(m):
                    if j == R[k][list(R[k]).index(i)+1]:
                        deltaT = 1 / C[k]
                        if isturn0:
                            deltaT *= 1.0 / (turn-turn0)
                        T[i, j] += deltaT
        #更新全局最优值
        local_best = min(C) #局部最优值
        bestk = list(C).index(local_best) #局部最优值对应的蚂蚁编号
        if local_best < global_best:
            global_best = local_best
            global_bestR = R[bestk].copy()
            turn0times = 0
            if isturn0:
                isturn1 = False

        else:
            turn0times += 1

        if turn0times == 5:
            turn0 = turn
            isturn0 = True

        turn += 1

        if turn1times == x and isturn1 and turn1num <= u:
            turn1num += 1
            sum_edge = 0
            num_edge = 0
            for i in range(n):
                for j in range(n):
                    if i == j or net[i, j] == float('inf'):
                        continue
                    sum_edge += 1.0 / net[i, j]
                    num_edge += 1;
            for i in range(n):
                for j in range(n):
                    if i == j or net[i, j] == float('inf'):
                        continue
                    T[i, j] = T[i, j] * num_edge * t0 * 1 / net[i, j] / sum_edge
            C = np.ones(m)    
            global_best = float('inf') 
            global_bestR = None 
            turn0times = 0 
            turn1times = 0  
            turn0 = None 
            isturn0 = False
            isturn1 = True 
            turn = 0           

    # '''输出'''
    # print('最短路径为：', global_best)
    # print('最短路径长度：', global_bestR)
    return global_best


def TanXin(n, net):
    R = np.ones(n+1, dtype=int)
    while True:
        for p in range(n+1):
            if p == 0:
                R[0] = ra.randint(0, n-1)
            elif p == n:
                R[n] = R[0]
            else:
                i = R[p-1]
                js = neighbor(n, net, i)
                js = list(set(js)-set(R[:p]))
                ljs = [net[i, j] for j in js]
                j = js[ljs.index(min(ljs))]
                R[p] = j
        C = sum([net[R[i], R[i+1]] for i in range(n)])
        if C != float('inf'):
            break
    return C

def neighbor(n, net, i):
    return [j for j in range(n) if j != i if net[i][j] != float('inf')]

if __name__ == '__main__':
    net = np.array([[float('inf'), 3, 1, 2],
                    [3, float('inf'), 5, 4],
                    [1, 5, float('inf'), 2],
                    [2, 4, 2, float('inf')]])
    TSP(net)

思路

TSP问题描述
TSP问题又称旅行商问题、货郎担问题，是数学领域的著名问题，也是计算机中的NP难问题。TSP问题的具体描述为：给定一系列城市和每对城市之间的距离，求解访问每一座城市一次并回到起始城市的最短回路。

蚁群算法的描述
蚁群算法的思想来自于自然界蚂蚁觅食，蚂蚁在寻找食物源时，会在路径上留下蚂蚁独有的路径标识——信息素，一条路径越短，该路径的信息素浓度就越高。蚂蚁会感知其他蚂蚁在各条路径上留下的信息素，更倾向于选择信息素浓度高的路径。随着时间的推移，越来越多的蚂蚁选择某一条路径，该路径的信息素浓度就越高，越能吸引蚂蚁，形成正反馈机制。最终大部分蚂蚁集中于最短的路径上。

基线模型
对于TSP问题，设蚁群中蚂蚁数量为(m)，结点数为(n)，结点i和结点j的距离为(net_{ij}), t时刻结点i和j之间连接的路径的信息素浓度为(T_{ij}(t))。初始时刻，所有边的信息素浓度设置为一个固定值，蚂蚁随机从不同的结点出发，按照一定概率选择路径。不妨设(P_{ij}^k(t))为在t时刻处于结点i的蚂蚁k选择结点j作为下一个结点的概率。蚂蚁选择下一个结点的概率收到两方面的影响，一是当前信息素浓度，二是访问下一个结点的长度，所以定义：

[P_{ij}^k(t)=frac{{[T_{ij}(t)]^a}cdot{[frac{1}{net_{ij}}]^b}}{sum_{jin allow_k}({[T_{ij}(t)]^a}cdot{[frac{1}{net_{ij}}]^b})} ]

其中(allow_k)为蚂蚁k下一步可以到达的结点集合，随着时间的推移，(allow_k)中的元素会越来越少，最终变为空集;(a)和(b)为重要程度因子
随着时间的推移，路径上信息素浓度会不断挥发，而蚂蚁行走也会增加路径上的信息素浓度，所以定义信息素的状态转移方程为：

[left{ egin{aligned} & T_{ij}(t+1) = (1-h) * T_{ij}(t) + sum_{k=1}^m(igtriangleup T_{ij}^k(t)) \ & igtriangleup T_{ij}^k(t) = 1 / C^k(t) end{aligned} ight. ]

其中(h)为信息素挥发率;(igtriangleup T_{ij}^k(t))为t时刻第k只蚂蚁在边(i, j)上增加的信息素浓度;(C^k(t))为在t时刻第k只蚂蚁经过路径的总长度。

改进模型
初始信息素浓度的改进：
在基线模型中，初始信息素浓度赋予相同的指，我认为这并不是一个好的方式，原因如下：
(1)一开始的信息素浓度不具有指导性，使得收敛速度下降
(2)可能使得蚂蚁一开始向着错误的路径前进，信息素在错误的路径上转移，容易陷入局部最优
所以我们根据路径大小给信息素浓度赋予不同的初值：

[T_{ij}(0) = Q cdot numEdge cdot frac{1/net_{ij}}{sum_{i,j}(1/net_{ij})} ]

其中，Q为初始浓度值，根据蚂蚁数量m除以贪心选择的路径长度得到;numEdge为网络的边数。该信息素初始化方法使得初始信息素浓度的平均值为Q,但是边的长度会影响具体的取值——边越短，初始信息素浓度越大。

信息素更新策略的改进：
在基线模型中，信息素浓度的更新策略是不变的。我们知道，运行一段时间后，蚁群会集中于某条路径，当前最优解往往就固定不变了，容易陷入局部最优。为了增加算法的全局搜索能力，减小局部最优发生的概率，需要一种更合理的信息更新策略。当迭代次数较小时，需要保证蚂蚁在路径上留下的信息素浓度足够高，以确保算法可以更快的收敛，而当算法收敛到一定程度时，如果多次迭代中算法所求得最有路径长度一直未发生变化，那么算法既有可能得到正确的结果，也有可能陷入局部最优，为了避免算法陷入局部最优的概率，需要适当增加信息素挥发速率和减少蚂蚁走过路径留下的信息素，避免信息素过于集中，使得蚂蚁能有更大概率去探索潜在的最优路径，提高算法的全局搜索能力。为了方便算法的描述，给出如下定义：
定义：如果存在某个时间节点t0,在一定固定的时长内，当前最优路径未发生变化，则称时间节点t0为稳定时刻点。

[left{ egin{aligned} &T_{ij}(t+1) = (1-h) * T_{ij}(t) + sum_{k=1}^m(igtriangleup T_{ij}^k(t)) & t < t0\ &T_{ij}(t+1) = frac{(1-h) * T_{ij}(t)}{t-t0} + frac{sum_{k=1}^m(igtriangleup T_{ij}^k(t))}{t-t0} & t >= t0 end{aligned} ight. ]

信息素矩阵的重置
为了进一步增加全家搜索的能力，在这里我们引入融断机制。当进入稳定时刻点t0后的一段时间内，全局最优解仍然没有发生变化，那么我们会将信息素矩阵重新初始化.

算法的流程图

算法评估
下面通过TSPLIB提供的三个数据集，来对比基线算法和改进算法的优劣。
数据集为：rand75.tsp eil76.tsp rat99.tsp
初始参数为: a=1, b=2, h=0.1 m=3
重复次数为：10

基线模型和改进模型的平均相对误差

数据集	基线模型	改进模型
数据集rand75.tsp	129%	57%
数据集eil76.tsp	180%	79%
数据集rat99.tsp	165%	73%

平均误差计算公式：(error = frac{predicted_value - measured_value}{measured_value})

改进模型较基线模型相对误差减少量

数据集	平均相对误差减少量
数据集rand75.tsp	72%
数据集eil76.tsp	101%
数据集rat99.tsp	92%

注意：数据的准确性不行，因为数据集大小不够，重复次数也不够，但是电脑不给力，时间不等人，只能这样了。

结论
可以看出，改进的模型相较与基线模型有一定的提升，但是依然存在一些问题：
(1)算法的时间消耗较大
(2)算法的结果相较于真实值还存在较大的相对误差。
(3)算法的参数没有缺少更多的实验
(4)具体的改进函数缺少更多的实验

参考

改进蚁群算法在TSP问题中的应用