蚁群算法简介及应用
蚂蚁几乎没有视力，但他们却能够在黑暗的世界中找到食物，而且能够找到一条从洞穴到食物的最短路径。它们是如何做到的呢？
简介
由来
蚁群算法是一种用来寻找优化路径的概率型算法。它由Marco Dorigo于1992年在他的博士论文中提出，其灵感来源于蚂蚁在寻找食物过程中发现路径的行为。
这种算法具有分布计算、信息正反馈和启发式搜索的特征，本质上是进化算法中的一种启发式全局优化算法。
思想
将蚁群算法应用于解决优化问题的基本思路为：用蚂蚁的行走路径表示待优化问题的可行解，整个蚂蚁群体的所有路径构成待优化问题的解空间。路径较短的蚂蚁释放的信息素量较多，随着时间的推进，较短的路径上累积的信息素浓度逐渐增高，选择该路径的蚂蚁个数也愈来愈多。最终，整个蚂蚁会在正反馈的作用下集中到最佳的路径上，此时对应的便是待优化问题的最优解。
应用
蚁群算法应用广泛，如旅行商问题(traveling salesman problem,简称TSP)、指派问题、Job-shop调度问题、车辆路径问题（vehicle routing problem）、图着色问题(graph coloring problem)和网络路由问题（network routing problem）等等，下面以TSP的求解为例。
scikit-opt库实现
import numpy as np from scipy import spatial import pandas as pd import matplotlib.pyplot as plt num_points = 25 points_coordinate = np.random.rand(num_points, 2) # generate coordinate of points distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean') def cal_total_distance(routine): num_points, = routine.shape return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)]) from sko.ACA import ACA_TSP aca = ACA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=50, max_iter=200, distance_matrix=distance_matrix) best_x, best_y = aca.run() # %% Plot fig, ax = plt.subplots(1, 2) best_points_ = np.concatenate([best_x, [best_x[0]]]) best_points_coordinate = points_coordinate[best_points_, :] ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r') pd.DataFrame(aca.y_best_history).cummin().plot(ax=ax[1]) plt.show()
Python实现
# -*- coding: utf-8 -*- import random import copy import time import sys import math import tkinter #//GUI模块 import threading from functools import reduce # 参数 ''' ALPHA:信息启发因子，值越大，则蚂蚁选择之前走过的路径可能性就越大 ，值越小，则蚁群搜索范围就会减少，容易陷入局部最优 BETA:Beta值越大，蚁群越就容易选择局部较短路径，这时算法收敛速度会 加快，但是随机性不高，容易得到局部的相对最优 ''' (ALPHA, BETA, RHO, Q) = (1.0,2.0,0.5,100.0) # 城市数，蚁群 (city_num, ant_num) = (50,50) distance_x = [ 178,272,176,171,650,499,267,703,408,437,491,74,532, 416,626,42,271,359,163,508,229,576,147,560,35,714, 757,517,64,314,675,690,391,628,87,240,705,699,258, 428,614,36,360,482,666,597,209,201,492,294] distance_y = [ 170,395,198,151,242,556,57,401,305,421,267,105,525, 381,244,330,395,169,141,380,153,442,528,329,232,48, 498,265,343,120,165,50,433,63,491,275,348,222,288, 490,213,524,244,114,104,552,70,425,227,331] #城市距离和信息素 distance_graph = [ [0.0 for col in range(city_num)] for raw in range(city_num)] pheromone_graph = [ [1.0 for col in range(city_num)] for raw in range(city_num)] #----------- 蚂蚁 ----------- class Ant(object): # 初始化 def __init__(self,ID): self.ID = ID # ID self.__clean_data() # 随机初始化出生点 # 初始数据 def __clean_data(self): self.path = [] # 当前蚂蚁的路径 self.total_distance = 0.0 # 当前路径的总距离 self.move_count = 0 # 移动次数 self.current_city = -1 # 当前停留的城市 self.open_table_city = [True for i in range(city_num)] # 探索城市的状态 city_index = random.randint(0,city_num-1) # 随机初始出生点 self.current_city = city_index self.path.append(city_index) self.open_table_city[city_index] = False self.move_count = 1 # 选择下一个城市 def __choice_next_city(self): next_city = -1 select_citys_prob = [0.0 for i in range(city_num)] #存储去下个城市的概率 total_prob = 0.0 # 获取去下一个城市的概率 for i in range(city_num): if self.open_table_city[i]: try : # 计算概率：与信息素浓度成正比，与距离成反比 select_citys_prob[i] = pow(pheromone_graph[self.current_city][i], ALPHA) * pow((1.0/distance_graph[self.current_city][i]), BETA) total_prob += select_citys_prob[i] except ZeroDivisionError as e: print ('Ant ID: {ID}, current city: {current}, target city: {target}'.format(ID = self.ID, current = self.current_city, target = i)) sys.exit(1) # 轮盘选择城市 if total_prob > 0.0: # 产生一个随机概率,0.0-total_prob temp_prob = random.uniform(0.0, total_prob) for i in range(city_num): if self.open_table_city[i]: # 轮次相减 temp_prob -= select_citys_prob[i] if temp_prob < 0.0: next_city = i break # 未从概率产生，顺序选择一个未访问城市 # if next_city == -1: # for i in range(city_num): # if self.open_table_city[i]: # next_city = i # break if (next_city == -1): next_city = random.randint(0, city_num - 1) while ((self.open_table_city[next_city]) == False): # if==False,说明已经遍历过了 next_city = random.randint(0, city_num - 1) # 返回下一个城市序号 return next_city # 计算路径总距离 def __cal_total_distance(self): temp_distance = 0.0 for i in range(1, city_num): start, end = self.path[i], self.path[i-1] temp_distance += distance_graph[start][end] # 回路 end = self.path[0] temp_distance += distance_graph[start][end] self.total_distance = temp_distance # 移动操作 def __move(self, next_city): self.path.append(next_city) self.open_table_city[next_city] = False self.total_distance += distance_graph[self.current_city][next_city] self.current_city = next_city self.move_count += 1 # 搜索路径 def search_path(self): # 初始化数据 self.__clean_data() # 搜素路径，遍历完所有城市为止 while self.move_count < city_num: # 移动到下一个城市 next_city = self.__choice_next_city() self.__move(next_city) # 计算路径总长度 self.__cal_total_distance() #----------- TSP问题 ----------- class TSP(object): def __init__(self, root, width = 800, height = 600, n = city_num): # 创建画布 self.root = root self.width = width self.height = height # 城市数目初始化为city_num self.n = n # tkinter.Canvas self.canvas = tkinter.Canvas( root, width = self.width, height = self.height, bg = "#EBEBEB", # 背景白色 xscrollincrement = 1, yscrollincrement = 1 ) self.canvas.pack(expand = tkinter.YES, fill = tkinter.BOTH) self.title("TSP蚁群算法(n:初始化 e:开始搜索 s:停止搜索 q:退出程序)") self.__r = 5 self.__lock = threading.RLock() # 线程锁 self.__bindEvents() self.new() # 计算城市之间的距离 for i in range(city_num): for j in range(city_num): temp_distance = pow((distance_x[i] - distance_x[j]), 2) + pow((distance_y[i] - distance_y[j]), 2) temp_distance = pow(temp_distance, 0.5) distance_graph[i][j] =float(int(temp_distance + 0.5)) # 按键响应程序 def __bindEvents(self): self.root.bind("q", self.quite) # 退出程序 self.root.bind("n", self.new) # 初始化 self.root.bind("e", self.search_path) # 开始搜索 self.root.bind("s", self.stop) # 停止搜索 # 更改标题 def title(self, s): self.root.title(s) # 初始化 def new(self, evt = None): # 停止线程 self.__lock.acquire() self.__running = False self.__lock.release() self.clear() # 清除信息 self.nodes = [] # 节点坐标 self.nodes2 = [] # 节点对象 # 初始化城市节点 for i in range(len(distance_x)): # 在画布上随机初始坐标 x = distance_x[i] y = distance_y[i] self.nodes.append((x, y)) # 生成节点椭圆，半径为self.__r node = self.canvas.create_oval(x - self.__r, y - self.__r, x + self.__r, y + self.__r, fill = "#ff0000", # 填充红色 outline = "#000000", # 轮廓白色 tags = "node", ) self.nodes2.append(node) # 显示坐标 self.canvas.create_text(x,y-10, # 使用create_text方法在坐标（302，77）处绘制文字 text = '('+str(x)+','+str(y)+')', # 所绘制文字的内容 fill = 'black' # 所绘制文字的颜色为灰色 ) # 顺序连接城市 #self.line(range(city_num)) # 初始城市之间的距离和信息素 for i in range(city_num): for j in range(city_num): pheromone_graph[i][j] = 1.0 self.ants = [Ant(ID) for ID in range(ant_num)] # 初始蚁群 self.best_ant = Ant(-1) # 初始最优解 self.best_ant.total_distance = 1 << 31 # 初始最大距离 self.iter = 1 # 初始化迭代次数 # 将节点按order顺序连线 def line(self, order): # 删除原线 self.canvas.delete("line") def line2(i1, i2): p1, p2 = self.nodes[i1], self.nodes[i2] self.canvas.create_line(p1, p2, fill = "#000000", tags = "line") return i2 # order[-1]为初始值 reduce(line2, order, order[-1]) # 清除画布 def clear(self): for item in self.canvas.find_all(): self.canvas.delete(item) # 退出程序 def quite(self, evt): self.__lock.acquire() self.__running = False self.__lock.release() self.root.destroy() print (u"\n程序已退出...") sys.exit() # 停止搜索 def stop(self, evt): self.__lock.acquire() self.__running = False self.__lock.release() # 开始搜索 def search_path(self, evt = None): # 开启线程 self.__lock.acquire() self.__running = True self.__lock.release() while self.__running: # 遍历每一只蚂蚁 for ant in self.ants: # 搜索一条路径 ant.search_path() # 与当前最优蚂蚁比较 if ant.total_distance < self.best_ant.total_distance: # 更新最优解 self.best_ant = copy.deepcopy(ant) # 更新信息素 self.__update_pheromone_gragh() print (u"迭代次数：",self.iter,u"最佳路径总距离：",int(self.best_ant.total_distance)) # 连线 self.line(self.best_ant.path) # 设置标题 self.title("TSP蚁群算法(n:随机初始 e:开始搜索 s:停止搜索 q:退出程序) 迭代次数: %d" % self.iter) # 更新画布 self.canvas.update() self.iter += 1 # 更新信息素 def __update_pheromone_gragh(self): # 获取每只蚂蚁在其路径上留下的信息素 temp_pheromone = [[0.0 for col in range(city_num)] for raw in range(city_num)] for ant in self.ants: for i in range(1,city_num): start, end = ant.path[i-1], ant.path[i] # 在路径上的每两个相邻城市间留下信息素，与路径总距离反比 temp_pheromone[start][end] += Q / ant.total_distance temp_pheromone[end][start] = temp_pheromone[start][end] # 更新所有城市之间的信息素，旧信息素衰减加上新迭代信息素 for i in range(city_num): for j in range(city_num): pheromone_graph[i][j] = pheromone_graph[i][j] * RHO + temp_pheromone[i][j] # 主循环 def mainloop(self): self.root.mainloop() #----------- 程序的入口处 ----------- if __name__ == '__main__': print (u""" -------------------------------------------------------- 程序：蚁群算法解决TPS问题程序 作者：许彬 日期：2015-12-10 语言：Python 2.7 -------------------------------------------------------- """) TSP(tkinter.Tk()).mainloop()
C实现
#include<bits/stdc++.h> using namespace std; const int INF = 0x3f3f3f3f; #define sqr(x) ((x)*(x)) #define eps 1e-8 string file_name; int type;// type == 1 全矩阵, type == 2 二维欧拉距离 int N;//城市数量 double **dis;//城市间距离 double **pheromone;//信息素 double **herustic;//启发式值 double **info;// info = pheromone ^ delta * herustic ^ beta double pheromone_0;//pheromone初始值，这里是1 / (avg * N)其中avg为图网中所有边边权的平均数。 int m;//种群数量 int delta, beta;//参数 double alpha; int *r1, *s, *r;//agent k的出发城市，下一个点，当前点。 int MAX, iteration;//最大迭代次数，迭代计数变量 set<int> empty, *J; struct vertex{ double x, y;// 城市坐标 int id;// 城市编号 int input(FILE *fp){ return fscanf(fp, "%d %lf %lf", &id, &x, &y); } }*node; typedef pair<int, int> pair_int; struct Tour{//路径 vector<pair_int> path;//path[i]，存储一条边(r,s) double L; void clean(){ L = INF; path.clear(); path.shrink_to_fit(); }//清空 void calc(){ L = 0; int sz = path.size(); for (int i = 0; i < sz; i ++){ L += dis[path[i].first][path[i].second]; } }//计算长度 void push_back(int x, int y){ path.push_back(make_pair(x, y)); } int size(){ return (int)path.size(); } int r(int i){ return path[i].first; } int s(int i){ return path[i].second; } void print(FILE *fp){ int sz = path.size(); for (int i = 0; i < sz; i ++){ fprintf(fp, "%d->", path[i].first + 1); } fprintf(fp, "%d\n", path[sz - 1].second + 1); } bool operator <(const Tour &a)const{ return L < a.L; }//重载 } *tour, best_so_far; double EUC_2D(const vertex &a, const vertex &b){ return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y)); } void io(){//输入 printf("input file_name and data type\n"); cin >> file_name >> type; FILE *fp = fopen(file_name.c_str(), "r"); fscanf(fp, "%d", &N); node = new vertex[N + 5]; dis = new double*[N + 5]; double tmp = 0; int cnt = 0; if (type == 1){ for (int i = 0; i < N; i ++){ dis[i] = new double[N]; for (int j = 0; j < N; j ++){ fscanf(fp, "%lf", &dis[i][j]); tmp += i != j ? dis[i][j] : 0;// i == j的时候dis不存在，所以不考虑。 cnt += i != j ? 1 : 0;// i == j的时候dis不存在，所以不考虑。 } } }else{ for (int i = 0; i < N; i ++) node[i].input(fp); for (int i = 0; i < N; i ++){ dis[i] = new double[N]; for (int j = 0; j < N; j ++){ dis[i][j] = EUC_2D(node[i], node[j]);// 计算距离 tmp += i != j ? dis[i][j] : 0;// i == j的时候 dis不存在，所以不考虑。 cnt += i != j ? 1 : 0;// i == j的时候dis不存在，所以不考虑。 } } } pheromone_0 = (double)cnt / (tmp * N);//pheromone初始值，这里是1 / (avg * N)其中avg为图网中所有边边权的平均数。 fclose(fp); return; } void init(){//初始化 alpha = 0.1;//evaporation parameter，挥发参数，每次信息素要挥发的量 delta = 1; beta = 6;// delta 和 beta分别表示pheromones和herustic的比重 m = N; pheromone = new double*[N + 5]; herustic = new double*[N + 5]; info = new double*[N + 5]; r1 = new int[N + 5]; r = new int[N + 5]; s = new int[N + 5]; J = new set<int>[N + 5]; empty.clear(); for (int i = 0; i < N; i ++){ pheromone[i] = new double[N + 5]; herustic[i] = new double[N + 5]; info[i] = new double[N + 5]; empty.insert(i); for (int j = 0; j < N; j ++){ pheromone[i][j] = pheromone_0; herustic[i][j] = 1 / (dis[i][j] + eps);//加一个小数eps，防止被零除 } } best_so_far.clean(); iteration = 0; MAX = N * N; } double power(double x, int y){//快速幂，计算x ^ y，时间复杂度O(logn),感兴趣可以百度 double ans = 1; while (y){ if (y & 1) ans *= x; x *= x; y >>= 1; } return ans; } void reset(){ tour = new Tour[m + 5]; for (int i = 0; i < N; i ++){ tour[i].clean(); r1[i] = i;//初始化出发城市， J[i] = empty; J[i].erase(r1[i]);//初始化agent i需要访问的城 市 r[i] = r1[i];//当前在出发点 } for (int i = 0; i < N; i ++) for (int j = 0; j < N; j ++){ info[i][j] = power(pheromone[i][j], delta) * power(herustic[i][j], beta); }//选择公式 } int select_next(int k){ if (J[k].empty()) return r1[k]; //如果J是空的，那么返回出发点城市 double rnd = (double)(rand()) / (double)RAND_MAX;//产生0..1的随机数 set<int>::iterator it = J[k].begin(); double sum_prob = 0, sum = 0; for (; it != J[k].end(); it ++){ sum += info[r[k]][*it]; }//计算概率分布 rnd *= sum; it = J[k].begin(); for (; it != J[k].end(); it ++){ sum_prob += info[r[k]][*it]; if (sum_prob >= rnd){ return *it; } }//依照概率选取下一步城市 } void construct_solution(){ for (int i = 0; i < N; i ++){ for (int k = 0; k < m; k ++){ int next = select_next(k);//选择下一步的最优决策 J[k].erase(next); s[k] = next; tour[k].push_back(r[k], s[k]); r[k] = s[k]; } } } void update_pheromone(){ Tour now_best; now_best.clean();//初始化 for (int i = 0; i < m; i ++){ tour[i].calc(); if (tour[i] < now_best) now_best = tour[i];//寻找当前迭代最优解 } if (now_best < best_so_far){ best_so_far = now_best;//更新最优解 } for (int i = 0; i < N; i ++) for (int j = 0; j < N; j ++) pheromone[i][j] *= (1 - alpha);//信息素挥发 int sz = now_best.size(); for (int i = 0; i < sz; i ++){ pheromone[now_best.r(i)][now_best.s(i)] += 1. / (double)now_best.L; pheromone[now_best.s(i)][now_best.r(i)] = pheromone[now_best.r(i)][now_best.s(i)];// 对称 }//更新信息素含量 } int main(){ srand((unsigned) time(0));//初始化随机种子 io(); init(); double last = INF; int bad_times = 0; for (; iteration < MAX; iteration ++){ if (bad_times > N) break;//进入局部最优 reset();//初始化agent信息 construct_solution();//对于所有的agent构造一个完整的tour update_pheromone();//更新信息素 printf("iteration %d:best_so_far = %.2f\n", iteration, best_so_far.L); if (last > best_so_far.L) last = best_so_far.L, bad_times = 0; else bad_times ++;//记录当前未更新代数，若迭代多次未更新，认为进入局部最优 } printf("best_so_far = %.2f\n", best_so_far.L);// 输出目标值 best_so_far.print(stdout);//输出路径 } 算例演示 例一 满秩矩阵式（type = 1） 输入文件格式为： File_name File_type salesman.in 1 5 0 1 2 2 3 2 0 3 4 2 3 2 0 4 1 3 4 5 0 5 2 4 1 4 0 输出结果为： opt_solution: 11 例二 二维坐标式（type = 2） 输入文件格式为： File_name File_type KroA100.tsp 2 100 1 1380 939 2 2848 96 3 3510 1671 4 457 334 5 3888 666 6 984 965 7 2721 1482 8 1286 525 9 2716 1432 10 738 1325 11 1251 1832 12 2728 1698 13 3815 169 14 3683 1533 15 1247 1945 16 123 862 17 1234 1946 18 252 1240 19 611 673 20 2576 1676 21 928 1700 22 53 857 23 1807 1711 24 274 1420 25 2574 946 26 178 24 27 2678 1825 28 1795 962 29 3384 1498 30 3520 1079 31 1256 61 32 1424 1728 33 3913 192 34 3085 1528 35 2573 1969 36 463 1670 37 3875 598 38 298 1513 39 3479 821 40 2542 236 41 3955 1743 42 1323 280 43 3447 1830 44 2936 337 45 1621 1830 46 3373 1646 47 1393 1368 48 3874 1318 49 938 955 50 3022 474 51 2482 1183 52 3854 923 53 376 825 54 2519 135 55 2945 1622 56 953 268 57 2628 1479 58 2097 981 59 890 1846 60 2139 1806 61 2421 1007 62 2290 1810 63 1115 1052 64 2588 302 65 327 265 66 241 341 67 1917 687 68 2991 792 69 2573 599 70 19 674 71 3911 1673 72 872 1559 73 2863 558 74 929 1766 75 839 620 76 3893 102 77 2178 1619 78 3822 899 79 378 1048 80 1178 100 81 2599 901 82 3416 143 83 2961 1605 84 611 1384 85 3113 885 86 2597 1830 87 2586 1286 88 161 906 89 1429 134 90 742 1025 91 1625 1651 92 1187 706 93 1787 1009 94 22 987 95 3640 43 96 3756 882 97 776 392 98 1724 1642 99 198 1810 100 3950 1558 输出结果为： best_known_solution: 21282
参考链接：
1. 百度百科-蚁群算法
4. CSDN-fanxin_i-蚁群算法原理及实现（Python）
5. 知乎-数据魔法师-十分钟快速get蚁群算法（附C代码）
个性签名：时间会解决一切
