蚁群算法(AntColonyOptimization,ACO)与TSP问题

对于小规模数据的TSP问题，我们可以使用动态规划快速的求解。

对于大规模数据的TSP问题，可以使用蚁群算法，模拟退火等近似算法进行求解。

蚁群算法是一种用来在图中寻找优化路径的机率型算法，最早Marco Dorigo提出。

它的灵感来源于蚁群寻找食物的过程，因为往往一只蚂蚁并没有太多“智能”的表现，而蚁群往往有“智能”的动作，比如大部分都趋向于食物

这是因为它们在移动的路径上会留下“信息素”，它们会更大概率沿着信息素更浓的路径行走，而路径越短，信息素就会越浓。

在解决TSP问题时，蚁群算法的流程如下：

1. 初始化算法
2. 派出一只蚂蚁，随机一个点作为蚂蚁的起点
3. 对于蚂蚁当前所在的点，根据所有后继边的信息素浓度以及长度根据公式计算出访问每个后继节点的概率
4. 用轮盘法随机出下一步走到哪个节点
5. 重复3,4步，直到所有节点都访问了一次
6. 计算回路长度，计算经过的边的信息素增加量
7. 所有边上的信息素都会蒸发一定比例使算法收敛，增加经过的边的信息素
8. 回到第2步，进行下一轮迭代

1.初始化问题：

初始化算法的时候需要设定一个初始信息素浓度，
需要注意这个浓度太大会导致新增的信息素没什么用，太小会导致过快结束算法，陷入局部最优解（类似于模拟退火的初始温度设定）
根据路径越长信息素浓度越低的规则，某只蚂蚁走过的回路长度为length的话，那么该蚂蚁对所走路径信息素浓度的贡献值为1/length
所以初始浓度设定为（M/length），M为蚂蚁数量

2.访问后继节点的概率计算：

算法的第三步需要计算当前节点到达后继节点的概率。

设定：α为信息素的重要程度，β为到达后继节点的边的长度，info[now][next]为信息素浓度，dist[now][next]为边长

那么到达一个后继节点的概率为： (info[now][next])^α/(dist[now][next])^β-----------公式一

假设蚂蚁现在处于A点，A点的后继节点为B，C，D

假设依据公式一所计算的结果，PB=0.5，PC=1.2，PD=0.2

那么：

选择节点B的概率为：pro_b = (PB)/(PB+PC+PD)=0.263

选择节点C的概率为：pro_c = (PC)/(PB+PC+PD)=0.632

选择节点D的概率为：pro_d = (PD)/(PB+PC+PD)=0.105

但是如果选择概率最高的节点作为下次抵达的节点，那么所有蚂蚁所作出的决策都是相同，这样就会失去探索新路径的机会，使得算法陷入停滞。

为了解决这个问题，我们采用轮盘法来做出抉择：

首先在[0,1]区间生成随机数RAND，然后进行迭代，RAND每次减去一个pro_*，当RAND减去一个pro_*后结果小于等于0，那么就选取节点*为下次抵达节点。

3.信号素的蒸发与增量

记录每只蚂蚁在跑图时，对于边edge[i][j]产生的信号素增量1/length[i][j]，。当蚁群中的所以蚂蚁跑完后，对于边edge[i][j]的信号素浓度info[i][j]按照一定比例进行蒸发操作

设蒸发系数为rho(rho一般取值为0.5~0.8，蒸发系数越大收敛越快，同时准确率降低)，那么修改后的info[i][j] = info[i][j]*(1-rho)

而后，将蚁群在跑图时，对于边edge[i][j]所产生的增量和phe[i][j]累加至info[i][j]

#include <bits/stdc++.h>
#include <iostream>
#define NS (25)
#define eps (1e-10)
#define M (80)
#define rho (0.5)  //信息素蒸发系数
#define alp (1)
#define bet (1)
#define Q (100)
#define Rand() ((double)rand() / RAND_MAX)

using namespace std;

int n, dis[NS][NS], path[NS], ans;

bool book[NS];

double info[NS][NS], phe[NS][NS], p[NS];

/*初始化过程：关于信息素的更新，
首先是初始化算法的时候需要设定一个初始信息素浓度，
这个浓度太大会导致新增的信息素没什么用，太小会导致过快结束算法，陷入局部最优解。
根据路径越长信息素浓度越低的规则，某只蚂蚁走过的回路长度为length的话，那么该蚂蚁对所走路径信息素浓度的贡献值为1/length
所以初始浓度设定为（Q*M/length）Q为常数用于精确计算结果，M为蚂蚁数量
*/
void init() //初始化算法
{
    int a = 1, len = 0;
    book[1] = 1;
    for (int c = 1; c < n; c += 1)
    {
        int mn = INT_MAX, nxt = 0;
        for (int i = 1; i <= n; i += 1)
            if (!book[i] && dis[a][i] < mn)
                mn = dis[a][i], nxt = i;
        len += dis[a][nxt], a = nxt, book[a] = 1;
    }
    len += dis[a][1], ans = len;
    for (int i = 1; i <= n; i += 1)
        for (int j = 1; j <= n; j += 1)
            info[i][j] = (double)Q * M / len;
}

inline double Pow(double a, int b)
{
    if(b==0)
        return 1;
    double res = 1;
    while(b > 0)
    {
        if(b & 1)
        {
            res = res*a;
        }
        a *= a;
        b /= 2;
    }
    return res;
}

/*
1.随机一个节点作为蚂蚁的起始节点
2.对于当前点，依据信息素浓度和与后继节点的边长计算访问该后继节点的概率
3.轮盘法随机出下一个到达的节点
4.重复2，3直到所有节点均被访问一次
*/
void run()
{
    int a = rand() % n + 1, s = a, len = 0;
    memset(book + 1, 0, sizeof(bool) * n), book[a] = 1;
    for (int c = 1; c < n; c += 1)
    {
        double tot = 0;
        for (int i = 1; i <= n; i += 1)
            if (book[i]) p[i] = 0;
            else
            {
                p[i] = Pow(info[a][i], alp) / Pow(dis[a][i], bet);
                tot += p[i];
            }
        if (tot < eps) return;
        for (int i = 1; i <= n; i += 1) p[i] /= tot;
        double r = Rand();
        for (int i = 1; i <= n; r -= p[i], i += 1)  //轮盘法
            if (!book[i] && r <= p[i])
            {
                len += dis[a][i], a = i, book[a] = 1, path[c] = i;
                break;
            }
    }
    // dt[][]为边上的信息素增量
    len += dis[a][s], phe[a][s] += (double)Q / len, ans = min(ans, len);
    for (int i = 1; i < n; i += 1)
        phe[s][path[i]] += (double)Q / len, s = path[i];
}

int main(int argc, char const* argv[])
{
    //cout << (0xfffffff) << endl;
    cin >> n;
    srand((unsigned) time(0));
    for (int i = 1; i <= n; i += 1)
        for (int j = 1; j <= n; j += 1)
            //IN(dis[i][j]);
            cin >> dis[i][j];
    init();
    for (int c = 1; c <= 700; c += 1)
    {
        for (int i = 1; i <= n; i += 1)
            for (int j = 1; j <= n; j += 1)
                phe[i][j] = 0;
        for (int j = 1; j <= M; j += 1) run();
        for (int i = 1; i <= n; i += 1)
            for (int j = 1; j <= n; j += 1)
                info[i][j] = info[i][j] * (1-rho) + phe[i][j];
    }
    printf("%d\n", ans);
    return 0;
}