找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
热搜: 活动 交友 discuz
查看: 24|回复: 0

蚁群算法 c++

[复制链接]

0

主题

0

回帖

26

积分

管理员

积分
26
发表于 2024-9-29 07:46:17 | 显示全部楼层 |阅读模式
什么是蚁群算法?

   蚁群系统(Ant System(AS)或Ant Colony System(ACS))是由意大利学者Dorigo、Maniezzo等人于20世纪90年代首先提出来的。他们在研究蚂蚁觅食的过程中,发现蚁群整体会体现一些智能的行为,例如蚁群可以在不同的环境下,寻找最短到达食物源的路径。

     后经进一步研究发现,这是因为蚂蚁会在其经过的路径上释放一种可以称之为“信息素(pheromone)”的物质,蚁群内的蚂蚁对“信息素”具有感知能力,它们会沿着“信息素”浓度较高路径行走,而每只路过的蚂蚁都会在路上留下“信息素”,这就形成一种类似正反馈的机制,这样经过一段时间后,整个蚁群就会沿着最短路径到达食物源了。


       由上述蚂蚁找食物模式演变来的算法,即是蚁群算法。这种算法具有分布计算、信息正反馈和启发式搜索的特征,本质上是进化算法中的一种启发式全局优化算法。

      最近几年,该算法在网络路由中的应用受到越来越多学者的关注,并提出了一些新的基于蚂蚁算法的路由算法。同传统的路由算法相比较,该算法在网络路由中具有信息分布式性、动态性、随机性和异步性等特点,而这些特点正好能满足网络路由的需要。


  蚁群算法应用广泛,如旅行商问题(traveling salesman problem,简称TSP)、指派问题、Job-shop调度问题、车辆路径问题(vehicle routing problem)、图着色问题(graph coloring problem)和网络路由问题(network routing problem)等等。


  1. #include<bits/stdc++.h>
  2. using namespace std;
  3. // const
  4. const int INF = 0x3f3f3f3f;
  5. #define sqr(x) ((x)*(x))
  6. #define eps 1e-8
  7. //variables
  8. string file_name;
  9. int type;// type == 1 全矩阵, type == 2 二维欧拉距离
  10. int N;//城市数量
  11. double **dis;//城市间距离
  12. double **pheromone;//信息素
  13. double **herustic;//启发式值
  14. double **info;// info = pheromone ^ delta * herustic ^ beta
  15. double pheromone_0;//pheromone初始值,这里是1 / (avg * N)其中avg为图网中所有边边权的平均数。
  16. int m;//种群数量
  17. int delta, beta;//参数
  18. double alpha;
  19. int *r1, *s, *r;//agent k的出发城市,下一个点,当前点。
  20. int MAX, iteration;//最大迭代次数,迭代计数变量
  21. set<int> empty, *J;
  22. struct vertex{
  23.                      double x, y;// 城市坐标
  24.                      int id;// 城市编号
  25.                      int input(FILE *fp){
  26.                      return fscanf(fp, "%d %lf %lf", &id,                          &x, &y);
  27. }
  28. }*node;
  29. typedef pair<int, int> pair_int;
  30. struct Tour{//路径
  31.        vector<pair_int> path;//path[i],存储一                     条边(r,s)
  32. double L;
  33.       void clean(){
  34.               L = INF;
  35.               path.clear();
  36.               path.shrink_to_fit();
  37.        }//清空
  38.       void calc(){
  39.               L = 0;
  40.               int sz = path.size();
  41.               for (int i = 0; i < sz; i ++){
  42.               L += dis[path[i].first][path[i].second];
  43.               }
  44.         }//计算长度
  45.       void push_back(int x, int y){
  46.               path.push_back(make_pair(x, y));
  47.         }
  48.       int size(){
  49.               return (int)path.size();
  50.         }
  51.       int r(int i){
  52.               return path[i].first;
  53.         }
  54.       int s(int i){
  55.               return path[i].second;
  56.         }
  57.       void print(FILE *fp){
  58.               int sz = path.size();
  59.               for (int i = 0; i < sz; i ++){
  60.               fprintf(fp, "%d->", path[i].first + 1);
  61.         }
  62.               fprintf(fp, "%d\n", path[sz - 1].second +                    1);
  63.         }
  64.       bool operator <(const Tour &a)const{
  65.                return L < a.L;
  66.         }//重载
  67. } *tour, best_so_far;
  68. double EUC_2D(const vertex &a, const vertex &b){
  69.       return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));
  70. }
  71. void io(){//输入
  72.       printf("input file_name and data type\n");
  73.       cin >> file_name >> type;
  74.       FILE *fp = fopen(file_name.c_str(), "r");
  75.       fscanf(fp, "%d", &N);
  76.       node = new vertex[N + 5];
  77.       dis = new double*[N + 5];
  78.       double tmp = 0;
  79.       int cnt = 0;
  80.       if (type == 1){
  81.                for (int i = 0; i < N; i ++){
  82.                     dis[i] = new double[N];
  83.                     for (int j = 0; j < N; j ++){
  84.                           fscanf(fp, "%lf", &dis[i][j]);
  85.                           tmp += i != j ? dis[i][j] : 0;// i == j的                           时候 dis不存在,所以不考虑。
  86.                           cnt += i != j ? 1 : 0;// i == j的时候                               dis不存在,所以不考虑。
  87.                     }
  88.                }
  89.         }else{
  90.               for (int i = 0; i < N; i ++)
  91.               node[i].input(fp);
  92.               for (int i = 0; i < N; i ++){
  93.                     dis[i] = new double[N];
  94.                     for (int j = 0; j < N; j ++){
  95.                           dis[i][j] = EUC_2D(node[i],                                         node[j]);// 计算距离
  96.                           tmp += i != j ? dis[i][j] : 0;// i == j的                           时候 dis不存在,所以不考虑。
  97.                           cnt += i != j ? 1 : 0;// i == j的时候                               dis不存在,所以不考虑。
  98.                      }
  99.               }
  100.        }
  101.         pheromone_0 =  (double)cnt / (tmp *                         N);//pheromone初始值,这里是1 / (avg * N)其中         avg为图网中所有边边权的平均数。
  102.         fclose(fp);
  103.         return;
  104. }
  105. void init(){//初始化
  106.         alpha = 0.1;//evaporation parameter,挥发参             数,每次信息素要挥发的量
  107.         delta = 1;
  108.         beta = 6;// delta 和 beta分别表示pheromones
  109.         和herustic的比重
  110.         m = N;
  111.         pheromone = new double*[N + 5];
  112.         herustic = new double*[N + 5];
  113.         info = new double*[N + 5];
  114.         r1 = new int[N + 5];
  115.         r = new int[N + 5];
  116.         s = new int[N + 5];
  117.         J = new set<int>[N + 5];
  118.         empty.clear();
  119.         for (int i = 0; i < N; i ++){
  120.                pheromone[i] = new double[N + 5];
  121.                herustic[i] = new double[N + 5];
  122.                info[i] = new double[N + 5];
  123.                empty.insert(i);
  124.                for (int j = 0; j < N; j ++){
  125.                       pheromone[i][j] = pheromone_0;
  126.                       herustic[i][j] = 1 / (dis[i][j] + eps);//加                         一个小数eps,防止被零除
  127.                 }
  128.          }
  129.          best_so_far.clean();
  130.          iteration = 0;
  131.          MAX = N * N;
  132. }
  133. double power(double x, int y){//快速幂,计算x ^ y,时间复杂度O(logn),感兴趣可以百度
  134.          double ans = 1;
  135.          while (y){
  136.                  if (y & 1) ans *= x;
  137.                  x *= x;
  138.                  y >>= 1;
  139.           }
  140.          return ans;
  141. }
  142. void reset(){
  143.          tour = new Tour[m + 5];
  144.          for (int i = 0; i < N; i ++){
  145.                 tour[i].clean();
  146.                 r1[i] = i;//初始化出发城市,
  147.                J[i] = empty;
  148.                J[i].erase(r1[i]);//初始化agent i需要访问的城                市
  149.                r[i] = r1[i];//当前在出发点
  150.          }
  151.          for (int i = 0; i < N; i ++)
  152.          for (int j = 0; j < N; j ++){
  153.               info[i][j] = power(pheromone[i][j], delta) *                 power(herustic[i][j], beta);
  154.          }//选择公式
  155. }
  156. int select_next(int k){
  157.          if (J[k].empty()) return r1[k]; //如果J是空的,那         么返回出发点城市
  158.         double rnd = (double)(rand()) /                                 (double)RAND_MAX;//产生0..1的随机数
  159.         set<int>::iterator it = J[k].begin();
  160.         double sum_prob = 0, sum = 0;
  161.         for (; it != J[k].end(); it ++){
  162.               sum += info[r[k]][*it];
  163.         }//计算概率分布
  164.         rnd *= sum;
  165.         it = J[k].begin();
  166.         for (; it != J[k].end(); it ++){
  167.              sum_prob += info[r[k]][*it];
  168.              if (sum_prob >= rnd){
  169.                     return *it;
  170.              }
  171.          }//依照概率选取下一步城市
  172. }
  173. void construct_solution(){
  174.        for (int i = 0; i < N; i ++){
  175.             for (int k = 0; k < m; k ++){
  176.                    int next = select_next(k);//选择下一步的                      最优决策
  177.                    J[k].erase(next);
  178.                    s[k] = next;
  179.                    tour[k].push_back(r[k], s[k]);
  180.                    r[k] = s[k];
  181.              }
  182.        }
  183. }
  184. void update_pheromone(){
  185.        Tour now_best;
  186.        now_best.clean();//初始化
  187.        for (int i = 0; i < m; i ++){
  188.             tour[i].calc();
  189.             if (tour[i] < now_best)
  190.             now_best = tour[i];//寻找当前迭代最优解
  191.        }
  192.        if (now_best < best_so_far){
  193.             best_so_far = now_best;//更新最优解
  194.        }
  195.        for (int i = 0; i < N; i ++)
  196.        for (int j = 0; j < N; j ++)
  197.        pheromone[i][j] *= (1 - alpha);//信息素挥发
  198.        int sz = now_best.size();
  199.        for (int i = 0; i < sz; i ++){
  200.             pheromone[now_best.r(i)][now_best.s(i)] +=             1. / (double)now_best.L;
  201.             pheromone[now_best.s(i)][now_best.r(i)] =               pheromone[now_best.r(i)][now_best.s(i)];//               对称
  202.        }//更新信息素含量
  203. }
  204. int main(){
  205.        srand((unsigned) time(0));//初始化随机种子
  206.        io();
  207.        init();
  208.        double last = INF;
  209.        int bad_times = 0;
  210.        for (; iteration < MAX; iteration ++){
  211.             if (bad_times > N) break;//进入局部最优
  212.             reset();//初始化agent信息
  213.             construct_solution();//对于所有的agent构造               一个完整的tour
  214.             update_pheromone();//更新信息素
  215.             printf("iteration %d:best_so_far = %.2lf\n",               iteration, best_so_far.L);
  216.             if (last > best_so_far.L)
  217.             last = best_so_far.L, bad_times = 0;
  218.             else bad_times ++;//记录当前未更新代数,若             迭代多次未更新,认为进入局部最优
  219.        }
  220.        printf("best_so_far = %.2lf\n", best_so_far.L);//        输出目标值
  221.        best_so_far.print(stdout);//输出路径
  222. }
复制代码
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Archiver|手机版|小黑屋|膜结构网

GMT+8, 2024-12-27 22:02 , Processed in 0.116208 second(s), 18 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表