|
什么是蚁群算法?
蚁群系统(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)等等。
- #include<bits/stdc++.h>
- using namespace std;
- // const
- const int INF = 0x3f3f3f3f;
- #define sqr(x) ((x)*(x))
- #define eps 1e-8
- //variables
- 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 = %.2lf\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 = %.2lf\n", best_so_far.L);// 输出目标值
- best_so_far.print(stdout);//输出路径
- }
复制代码 |
|