標籤：

【智能優化演算法?愛鑽研】-3十分鐘快速get蟻群演算法（附代碼）

05-10

之前分享了TSP的動態規劃解法，本期來介紹它的另一種解法——蟻群演算法。

內容提要：

＊什麼是蟻群演算法

＊蟻群演算法演練

＊演算法補充筆記

什麼是蟻群演算法？

蟻群系統(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）等等。

下面我們同之前推文一樣，以TSP的求解為例演練蟻群演算法的應用。

關於TSP問題，如果還有疑問，請參考之前的推文:

a href="mp.weixin.qq.com/s?"> 「乾貨|十分鐘教你用動態規劃演算法解Travelling

Salesman Problem（TSP）問題，附代碼……」。

關於求解TSP的蟻群演算法可以參考文章:

Ant colony

system: a cooperative learning approach to the traveling salesman problem,M.

Dorigo, L.M. Gambardella, IEEE Transactions on Evolutionary Computation,Volume:

1, Issue: 1, Apr 1997, pages 53 - 66。

蟻群演算法求解TSP

1. TSP建模

2. 蟻群演算法

附. 蟻群演算法相關代碼

先放上一波嚴肅的偽代碼分析：

接下來分享一波*真*代碼和算例 ↓

代碼（C++）

#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
", 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
");

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
",

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
", 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

注意：

使用本程序的時候只需要建立一個上述文件名的文檔，放在與源程序同目錄下面，並運行程序，輸入文件名以及數據類型。

例如：

運行程序之後會遇到以下提升：

input file_name and data type

只需要輸入

KroA100.tsp 2

即可得到一個啟發解以及相應路徑

更多的數據可以從TSPLIB下載。

上述代碼僅供分享交流學慣用，如有需要複製下面鏈接自取 ↓ ↓ ↓

paste.ubuntu.com/253532

演算法補充筆記

實際實驗中發現，當螞蟻在一條路徑上覓食很久時，放置一個近的食物基本沒有效果，這可以理解為當一隻螞蟻找到一條路徑時，在經過很長時間後大多數螞蟻都選擇了這條路徑，這時，突然有一隻螞蟻找到了較近的食物，但因為時間過得太久，兩條路徑上濃度相差太大（濃度越大，被選擇的概率就越大），整個系統基本已經停滯了，陷入了局部最優。所以簡單的螞蟻系統是存在一些問題的，如：

搜索到一定程度，會出現停滯狀態，陷入局部最優的情況

↓ ↓ ↓

經過小編的十分鐘·蟻群演算法·快速·真·攻略的分享，相信現在世界上又多了一批建設 和·諧·世·界 的演算法master！

如果大家對蟻群演算法 及 文中所敘內容 還有疑問或想要交流心得建議，歡迎在推文下留言溝通！

—end—

編輯：謝良楨（1922193128@qq.com）

賀興（hexing15@gmail.com）

代碼：賀興（hexing15@gmail.com）

指導老師：秦時明岳（professor.qin@qq.com）

如有疑問，歡迎諮詢~

更多最新數據分析相關原創文章，請關注微信公眾號：數據魔術師


推薦閱讀：

※二叉樹中相關遞歸求解(c,c++數據結構，演算法設計與分析)
※028 Implement strStr()[E]
※Q-learning 和 DQN
※MATLAB實用技巧&快捷方式&編程方法
※014 Longest Common Prefix[E]

TAG:演算法設計 | 蟻群演算法 |