C++蚁群算法，带图，超详细

admin

1. //main.cpp
2. #include"Head.h"
4. int main() {
5.         //更新随机种子
6.         srand((unsigned)time(NULL));
8.         Graph* graph = new Graph;
10.         graph->readCity();//启用文件读取，注销掉则为随机构造城市
12.         graph->innitial_Graph();
13.         AntColony* ant_colony = new AntColony;
14.         ant_colony->innitial_AntColony();
16.         for (int NC_now = 0; NC_now < NC_MAX; ++NC_now) {
17.                 for (int i = 1; i < N; ++i) {
18.                         ant_colony->ChooseNextCity(graph, i);
19.                 }
20.                 ant_colony->Update(graph,NC_now);
21.                 ant_colony->innitial_AntColony();
22.         }
23.        
24.         cout << "The BestLength=" << BestLength << endl;
26.         ShowByPython(graph);
28.         int last, next;
29.         double sum = 0;
30.         for (int i = 1; i < N; ++i) {
31.                 last = BestPath[i-1];
32.                 next = BestPath[i];
33.                 cout << "Travel " << i << " : " << last << "   \t-> " << next << " \tdistante: " << graph->Length[last][next] << endl;
34.                 sum += graph->Length[last][next];
35.         }
36.         last = BestPath[N - 1];
37.         next = BestPath[0];
38.         cout << "Travel " << N << " : " << last << "   \t-> " << next << " \tdistante: " << graph->Length[last][next] << endl;
39.         sum += graph->Length[last][next];
40.         cout << "最优解路程和 = " << sum << endl;
41.         cout << "TSPlib上的已知最优解 :" << graph->KnowBest << endl;
42.         cout << "与已知最优解的偏差为 :" << ((sum - graph->KnowBest) / graph->KnowBest) * 100 << "%" << endl;//最后两句必须启用文件读取才有效
44.         delete graph;
45.         delete ant_colony;
46.         return 0;
47. }

复制代码

1. //Head.h
2. #include<iostream>
3. #include<cstdlib>
4. #include<string>
5. #include<cmath>
6. #include<fstream>
7. #include<Python.h>
8. using namespace std;
9. #define PATH "E:\\TSPlib\\ch150.txt"//TSPlib数据文件的路径
11. #define N 150                //城市数量
12. #define M 100                //蚂蚁数量
13. #define NC_MAX 1200        //迭代次数
14. #define α 1                //信息素浓度重要性
15. #define β 1                //启发式重要性
16. #define Q 100                //蚂蚁所携带信息素
17. #define ρ 0.2                //蒸发系数
19. //注意：如果从文件中读取City，目前需要手动调整城市数量N，要与文件中保持一致
20. //并且修改 PATH 路径与启用main函数中的readCity函数
22. //以下数据结构用于最后分析
23. int BestPath_PerRound[NC_MAX][N] = { 0 };//记录每轮的最优路径
24. double BestLength_PerRound[NC_MAX] = { 0 };//记录每轮最优路径的长度
25. int BestPath[N] = { 0 };//全局最短路径（后一轮迭代未必一定比前一轮更优，只是总体是这样，不绝对），不保存每轮的信息
26. double BestLength = 1000000;//全局最短路径长度
27. double AverageLength_PerRound[NC_MAX] = { 0 };//记录每轮所有蚂蚁所走路径的平均长度
29. class City {
30. public:
31.         //无参构造函数，创建城市时默认在[0,999]内随机获得x,y坐标
32.         City() {
33.                 x = rand() % 1000;
34.                 y = rand() % 1000;
35.         }
36.         double x;
37.         double y;
38. };
39. class Graph {
40. public:
41.         City citys[N];//无参构造函数创建N个城市
42.         double Length[N][N] = { 0 };//路径长度表
43.         double tao[N][N] = { 0 };//信息素浓度表
44.         double Heuristic[N][N] = { 0 };//启发式表，将城市之间距离的倒数作为启发函数
45. public:
46.         //使用fstream读取文件中的城市坐标，并创建城市
47.         int KnowBest;//读取文件时存放已知最优解
48.         void readCity() {
49.                 ifstream CityFile;
50.                 CityFile.open(PATH, ios::in);
51.                 if (!CityFile.is_open()) {
52.                         cout << "Open file failure" << endl;
53.                         exit(0);
54.                 }
55.                 CityFile >> KnowBest;
56.                 cout << "The known best result is :" << KnowBest << endl;
57.                 int a;double b[2];
58.                 while (!CityFile.eof()) {
59.                         CityFile >> a >> b[0] >> b[1];
60.                         citys[a - 1].x = b[0];
61.                         citys[a - 1].y = b[1];
62.                 }
63.                 CityFile.close();
64.         }
65.         //初始化城市图，并计算各城市之间的距离
66.         void innitial_Graph() {
67.                 double x = 0, y = 0;
68.                 for (int i = 0; i < N; ++i) {
69.                         for (int j = 0; j < N; ++j) {
70.                                 x = citys[i].x - citys[j].x;
71.                                 y = citys[i].y - citys[j].y;
72.                                 Length[i][j] = sqrt(pow(x, 2) + pow(y, 2));//就是两点间距离公式
73.                                 if (i == j) {
74.                                         Length[i][j] = 10000;//默认城市到自己的距离为10000，防止启发式表出现无穷大，并不会计算城市自己到自己的概率
75.                                 }
76.                                 else {
77.                                         if (x == 0 && y == 0) {
78.                                                 //这种情况是城市坐标完全重合，如果遇到这种情况则不能继续计算了，否则i与j之间的距离为0，启发式就会为无穷大
79.                                                 //随之p的分子和分母都会变成无穷大，无穷大除无穷大...
80.                                                 cout << "innitial_Graph ERROR!!! City overlapping, please retry." << endl;
81.                                                 cout << "City :" << i << " and City :" << j << endl;
82.                                                 exit(0);
83.                                         }
84.                                 }
85.                                 Heuristic[i][j] = 1 / Length[i][j];//将距离的倒数作为启发式
86.                                 tao[i][j] = 1;//信息素tao不能初始化为0
87.                         }
88.                 }
89.         }
90. };
91. class AntColony {
92. public:
93.         int tabu[M][N] = { 0 };//禁忌列表，同时记录了蚂蚁的路径，每个元素取值为[0,N-1]
94.         int allow[M][N] = { 0 };//允许列表，完全依赖于禁忌列表，只是为了算法方便一些而设置，无法记录路径，每个元素取值为[0,1]
95.         double probability[M][N] = { 0 };//概率表，存放各蚂蚁选择下一作城市的概率
96.         double cumulative_probability[M][N] = { 0 };//累积概率表，用于轮盘赌算法随机选择城市
97. public:
99.         //每轮迭代前都应该被调用一次
100.         void innitial_AntColony() {
101.                 //初始化各列表
102.                 for (int i = 0; i < M; ++i) {
103.                         for (int j = 0; j < N; ++j) {
104.                                 tabu[i][j] = -1;
105.                                 allow[i][j] = 1;
106.                                 probability[i][j] = 0;
107.                                 cumulative_probability[i][j] = 0;
108.                         }
109.                 }
110.                 //为每一只蚂蚁随机一个初始城市，并将其加入禁忌列表
111.                 for (int i = 0; i < M; i++) {
112.                         tabu[i][0] = rand() % N;
113.                         allow[i][tabu[i][0]] = 0;
114.                 }
115.         }
117.         //求取概率表和累积概率表，然后根据轮盘赌算法选择下一个城市，每轮迭代应该被调用N-1次（初始化时已经有一个城市了）
118.         void ChooseNextCity(Graph* graph, int times) {//times表示此轮迭代已经是第几次调用这个函数了，times应该从1开始，到N-1次结束
120.                 //求概率表
121.                 double sum = 0;
122.                 int city_now = -1;
123.                 for (int i = 0; i < M; ++i) {
124.                         sum = 0;
125.                         city_now = tabu[i][times - 1];//蚂蚁i目前在城市city_now
126.                         for (int j = 0; j < N; ++j) {
127.                                 if (allow[i][j] == 1) {
128.                                         probability[i][j] = pow(graph->tao[city_now][j], α) * pow(graph->Heuristic[city_now][j], β);
129.                                         sum += probability[i][j];
130.                                 }
131.                                 else {
132.                                         probability[i][j] = 0;
133.                                         sum += 0;
134.                                 }
135.                         }
136.                         //上面求出的probability表并不是真正的概率表，而只是概率的分子，下面求出真正的概率表
137.                         for (int j = 0; j < N; ++j) {
138.                                 probability[i][j] = probability[i][j] / sum;
139.                         }
140.                 }
142.                 //用概率表求累积概率表
143.                 for (int i = 0; i < M; ++i) {
144.                         cumulative_probability[i][0] = probability[i][0];
145.                         for (int j = 1; j < N; ++j) {
146.                                 cumulative_probability[i][j] = cumulative_probability[i][j - 1] + probability[i][j];
147.                         }
148.                 }
150.                 //依据累积概率表，用轮盘赌算法为每只蚂蚁选择下一个城市
151.                 double p = 0;
152.                 for (int i = 0; i < M; ++i) {
153.                         p = rand() % 1000;
154.                         p /= 1000;
155.                         for (int j = 0; j < N; ++j) {
156.                                 if (p <= cumulative_probability[i][j]) {
157.                                         //如果满足此式，则让蚂蚁i前往城市j（下面更新tabu表的操作就是从逻辑上把蚂蚁移动到城市j）
158.                                         tabu[i][times] = j;
159.                                         allow[i][j] = 0;
160.                                         break;
161.                                 }
162.                         }
163.                 }
164.         }
166.         //更新函数每次迭代后都应该被调用一次，用于更新信息素表graph->tabu，并记录本轮迭代的相关信息（最优路径，最优路径的长度，所有蚂蚁所走路径的平均长度）
167.         //NC_now表示目前的迭代次数，应该从0开始，到NC_MAX-1结束，依次调用此函数
168.         void Update(Graph* graph, int NC_now) {
169.                 double delta_ktao = 0;        //用于保存当前蚂蚁留在所经边的信息素
170.                 double delta_tao[N][N] = { 0 };        //此轮迭代的信息素增量表，表示本轮信息素的增量，配合蒸发系数更新信息素表
171.                 double sum_Length[M] = { 0 };        //保存此轮每只蚂蚁所经过的路径长度
173.                 int last_city, next_city;
174.                 //遍历tabu表，计算每只蚂蚁的路径长度，存放在sum_Length[]中
175.                 for (int i = 0; i < M; ++i) {
176.                         for (int j = 1; j < N; ++j) {
177.                                 last_city = tabu[i][j - 1];
178.                                 next_city = tabu[i][j];
179.                                 sum_Length[i] += graph->Length[last_city][next_city];
180.                         }
181.                         //还要再加上终点到起点的距离
182.                         last_city = tabu[i][N - 1];
183.                         next_city = tabu[i][0];
184.                         sum_Length[i] += graph->Length[last_city][next_city];
185.                 }
186.                 //遍历tabu表，计算信息素增量表delta_tao[][]
187.                 for (int i = 0; i < M; ++i) {
188.                         delta_ktao = Q / sum_Length[i];
189.                         for (int j = 1; j < N; ++j) {
190.                                 last_city = tabu[i][j - 1];
191.                                 next_city = tabu[i][j];
192.                                 delta_tao[last_city][next_city] += delta_ktao;
193.                         }
194.                 }
195.                 //利用信息素增量表和蒸发系数更新信息素表tao[][]
196.                 for (int i = 0; i < N; ++i) {
197.                         for (int j = 0; j < N; ++j) {
198.                                 graph->tao[i][j] = graph->tao[i][j] * (1 - ρ) + delta_tao[i][j];
199.                         }
200.                 }
202.                 //计算本轮最优路径，最优路径的长度，所有蚂蚁所走路径的平均长度
203.                 int flag = 0;
204.                 double min = 1000000, sum = 0;        //min的初始值必须是一个足够大的值
205.                 for (int i = 0; i < M; ++i) {
206.                         sum += sum_Length[i];
207.                         if (min > sum_Length[i]) {
208.                                 min = sum_Length[i];
209.                                 flag = i;//标记最好的蚂蚁
210.                         }
211.                 }
212.                 //记录本轮信息至全局变量中
213.                 for (int i = 0; i < N; ++i) {
214.                         BestPath_PerRound[NC_now][i] = tabu[flag][i];
215.                 }
216.                 BestLength_PerRound[NC_now] = min;
217.                 AverageLength_PerRound[NC_now] = sum / M;
220.                 //更新全局最优路径和其长度
221.                 if (BestLength > min) {
222.                         for (int i = 0; i < N; ++i) {
223.                                 BestPath[i] = BestPath_PerRound[NC_now][i];
224.                         }
225.                         BestLength = min;
226.                 }
228.                 //用2-opt局部搜索优化本轮最优路径并用信息素强化
229.                 {
230.                         int TempPath[N];
231.                         int flag = false;
232.                         int t;
233.                         for (int i = 0; i < N; ++i) {
234.                                 TempPath[i] = BestPath_PerRound[NC_now][i];
235.                         }
236.                         for (int a = 0; a < N-1; ++a) {
237.                                 //如果路径被优化为新路径，则再对新路径从头来一次局部搜索
238.                                 if (flag == true) {
239.                                         a = 0;
240.                                         flag = false;
241.                                 }
242.                                 for (int b = a + 1; b < N; ++b) {
243.                                         //逆序a~b的路径
244.                                         for (int i = a, j = b; i < j; ++i, --j) {
245.                                                 t = TempPath[i];
246.                                                 TempPath[i] = TempPath[j];
247.                                                 TempPath[j] = t;
248.                                         }
249.                                         //重新求优化后的路径长度，加上终点到起点的距离
250.                                         sum = 0;
251.                                         for (int i = 1; i < N; ++i) {
252.                                                 last_city = TempPath[i - 1];
253.                                                 next_city = TempPath[i];
254.                                                 sum += graph->Length[last_city][next_city];
255.                                         }
256.                                         last_city = TempPath[N - 1];
257.                                         next_city = TempPath[0];
258.                                         sum += graph->Length[last_city][next_city];
259.                                         //如果此解更优则更新本轮最优解
260.                                         if (sum < BestLength_PerRound[NC_now]) {
261.                                                 flag = true;//如果路径被更新，则局部搜索从0开始
262.                                                 //先更新平均路径再更新最优路径
263.                                                 AverageLength_PerRound[NC_now] = (AverageLength_PerRound[NC_now] * M - BestLength_PerRound[NC_now] + sum) / M;
264.                                                 BestLength_PerRound[NC_now] = sum;
265.                                                 for (int i = 0; i < N; ++i) {
266.                                                         BestPath_PerRound[NC_now][i] = TempPath[i];
267.                                                 }
268.                                                 //如果比全局最优解还优，则更新全局最优解
269.                                                 if (sum < BestLength) {
270.                                                         for (int i = 0; i < N; ++i) {
271.                                                                 BestPath[i] = TempPath[i];
272.                                                         }
273.                                                         BestLength = sum;
274.                                                 }
275.                                         }
276.                                 }
277.                         }
278.                         //使用2opt局部搜索获得更强的解之后，给最优路径追加信息素奖励
279.                         double reward = Q / sum;
280.                         for (int i = 1; i < N; ++i) {
281.                                 last_city = BestPath[i-1];
282.                                 next_city = BestPath[i];
283.                                 graph->tao[last_city][next_city] += reward;
284.                         }
285.                 }
286.         }
287. };
289. //此函数的目的是为了用图展示出蚁群算法的结果，用C++调用python的matplotlib库
290. bool ShowByPython(Graph* graph) {
292.         //C++中初始化python环境
293.         Py_Initialize();
294.         if (!Py_IsInitialized()) {
295.                 std::cout << "Py_Initialize Error!" << std::endl;
296.                 return false;
297.         }
299.         try {
300.                 //执行python语句
301.                 PyRun_SimpleString("import sys");
302.                 PyRun_SimpleString("import matplotlib.pyplot as plt");
303.                 string importPy = "sys.argv = ['suibiantian']";
304.                 PyRun_SimpleString(importPy.c_str());
306.                 double x, y;
307.                 PyRun_SimpleString("plt.figure(num=1)");
308.                 PyRun_SimpleString("plt.title('Ant Colony algorithm for solving TSP')");//蚁群算法解决TSP问题
309.                 PyRun_SimpleString("plt.xlabel('x')");
310.                 PyRun_SimpleString("plt.ylabel('y')");
311.                 importPy= "plt.xlabel('Shortest=" + to_string(BestLength) + "')";
312.                 PyRun_SimpleString(importPy.c_str());
313.                 for (int i = 0; i < N; ++i) {
314.                         x = graph->citys[i].x;//Graph[i].get_x();
315.                         y = graph->citys[i].y;
316.                         importPy = "plt.scatter(" + to_string(x) + "," + to_string(y) + ")";
317.                         PyRun_SimpleString(importPy.c_str());
318.                 }
320.                 int last_city, next_city;
321.                 double X1, X2, Y1, Y2;
322.                 string str;
323.                 //把全局最优路径按顺序画出来,画在第一副图上
324.                 PyRun_SimpleString("plt.figure(num=1)");
325.                 for (int i = 1; i < N; ++i) {
326.                         last_city = BestPath[i - 1];
327.                         next_city = BestPath[i];
328.                         X1 = graph->citys[last_city].x;
329.                         Y1 = graph->citys[last_city].y;
330.                         X2 = graph->citys[next_city].x;
331.                         Y2 = graph->citys[next_city].y;
332.                         str = "plt.plot([" + to_string(X1) + "," + to_string(X2) + "],[" + to_string(Y1) + "," + to_string(Y2) + "])";
333.                         PyRun_SimpleString(str.c_str());
334.                 }
335.                 //把终点到起点的线也画上
336.                 last_city = BestPath[N - 1];
337.                 next_city = BestPath[0];
338.                 X1 = graph->citys[last_city].x;
339.                 Y1 = graph->citys[last_city].y;
340.                 X2 = graph->citys[next_city].x;
341.                 Y2 = graph->citys[next_city].y;
342.                 str = "plt.plot([" + to_string(X1) + "," + to_string(X2) + "],[" + to_string(Y1) + "," + to_string(Y2) + "])";
343.                 PyRun_SimpleString(str.c_str());
345.                 //下面展示随着迭代进行，每轮最优路径长度的变化，这是第二幅图
346.                 double last = 0, next = 0;
347.                 //string str;
348.                 PyRun_SimpleString("plt.figure(num=2)");
349.                 PyRun_SimpleString("plt.title('ShortestLength Convergence graph')");//最优路径收敛图
350.                 PyRun_SimpleString("plt.xlabel('times of iterations')");
351.                 PyRun_SimpleString("plt.ylabel('BestLength')");
352.                 for (int i = 1; i < NC_MAX; ++i) {
353.                         last = BestLength_PerRound[i - 1];
354.                         next = BestLength_PerRound[i];
355.                         str = "plt.plot([" + to_string(i - 1) + "," + to_string(i) + "],[" + to_string(last) + "," + to_string(next) + "])";
356.                         PyRun_SimpleString(str.c_str());
357.                 }
359.                 //下面展示随着迭代进行，每轮平均路径长度的变化，这是第三幅图
360.                 PyRun_SimpleString("plt.figure(num=3)");
361.                 PyRun_SimpleString("plt.title('AverageLength Convergence graph')");//平均路径收敛图
362.                 PyRun_SimpleString("plt.xlabel('times of iterations')");
363.                 PyRun_SimpleString("plt.ylabel('AverageLength')");
364.                 for (int i = 1; i < NC_MAX; ++i) {
365.                         last = AverageLength_PerRound[i - 1];
366.                         next = AverageLength_PerRound[i];
367.                         str = "plt.plot([" + to_string(i - 1) + "," + to_string(i) + "],[" + to_string(last) + "," + to_string(next) + "])";
368.                         PyRun_SimpleString(str.c_str());
369.                 }
371.                 //展示上面描绘的三幅图
372.                 PyRun_SimpleString("plt.show()");
373.         }
374.         catch (...) {
375.                 PyErr_Print();
376.                 PyErr_Clear();
377.                 Py_Finalize();
378.                 return false;
379.         }
380.         Py_Finalize();
381.         return true;
382. }

复制代码