嗯哼，時隔半年，再次有時間整理關於組合優化問題——旅行商問題（Traveling Salesman Problem, TSP），這次采用的是經典遺傳算法（Genetic Algorithm， GA）進行求解，利用C++語言進行編程實現。關於TSP問題以及GA的簡單介紹，可參見我的另一篇文章：Java版GA_TSP（我的第一個Java程序）。
各種啟發式算法的整體框架大致都由以下幾個操作組成：（1）初始解的產生；（2）解的評價（評價函數）；（3）擾動算子；此外，還可以加上程序原始數據的導入等操作。這些操作是多數啟發式算法所通用的算子，基為此，此次在采用C++進行實現的時候，采用一個通用的 HeuristicOperator.h 頭文件以及對應的 HeuristicOperator.cpp 類文件對這些操作進行集中放置，造好輪子，方便以后取用。
Ps:指針玩不太轉，只好繞道過去了~~~~
表1 HeuristicOperator中函數功能清單
| 編號 |
功能 |
記號 |
| 函數1 |
獲取客戶點坐標 |
getCoord() |
| 函數2 |
獲取距離矩陣 |
getDM() |
| 函數3 |
獲取初始解 |
getInitS() |
| 函數4 |
解的評價 |
Eval() |
| 函數5 |
搜索范圍內最優評價值及其相對應的位置 |
bestS() |
| 函數6 |
產生Sharking操作位置 |
RandPosition() |
| 函數7 |
交換算子（swap） |
Swap() |
| 函數8 |
翻轉算子（flip） |
Flip() |
| 函數9 |
插入算子（insert） |
Insert() |
| 注：函數5可以直接用 STL 中vector的操作函數max_element、min_element實現類似的功能，寫的時候一時沒有想起來，就自己造了個輪子。 |
||
之前在采用Matlab以及Java實現GA to solve TSP 時，考慮到程序的運行效率，通常會選用 array 來放置各種數據，眾所周知，array的特點就是固定分配的連續物理地址進行數據的存儲，然而對於不定長度的數據進行存儲時，一般的方式是采用“多次少量”，即預先分配一定內存空間，等不夠用了再分配一定的內存空間（多說一句，Matlab 中還可以采用以下兩種方式：用 array=[]; 以及用cell實現存儲不定長度數組，但效率不高）。而在C++ STL 中有一種神奇的數據類型vector容器，它既有着 array 的連續內存分配方式，又能不用指定數據存儲長度，對於一組不同規模的數據集進行測試時，再也不用擔心使用array時提醒必須預分配確定的存儲空間了~~
以下為HeuristicOperator.h頭文件：
1 #pragma once
2 #include<iostream>
3 #include<vector>
4 #include <algorithm> // std::shuffle
5 #include <random> // std::default_random_engine
6 #include<chrono>
7 using namespace std; 8
9 class HeuristicOperator { 10 public: 11 vector<vector<double>> getCoord(void); //函數1：獲取坐標函數
12 vector<vector<double>> getDM(vector<vector<double>> Coord); //函數2：獲取距離矩陣函數
13 vector<int> getInitS(int n); //函數3：獲取初始解函數
14 double Eval(vector<int> S, vector<vector<double>> DM, int n); //函數4：評價函數
15
16 vector<double> bestS(vector<double> Eval, int Length); //函數5：搜索范圍內最優評價值及其相應的位置函數
17
18 vector<int> RandPosition(int n); //函數6：產生Sharking操作位置函數
19 vector<int> Swap(vector<int> S, vector<int> RP); //函數7：交換算子
20 vector<int> Flip(vector<int> S, vector<int> RP); //函數8：翻轉算子
21 vector<int> Insert(vector<int> S, vector<int> RP); //函數9：插入算子
22 };
對應的HeuristicOperator.cpp類文件比較容易實現，在此不再贅述。本文所用算例為31城市的TSP問題，與Java版遺傳算法求解TSP求解算例一致，具體數據如下：
1 1304 2312
2 3639 1315
3 4177 2244
4 3712 1399
5 3488 1535
6 3326 1556
7 3238 1229
8 4196 1004
9 4312 790
10 4386 570
11 3007 1970
12 2562 1756
13 2788 1491
14 2381 1676
15 1332 695
16 3715 1678
17 3918 2179
18 4061 2370
19 3780 2212
20 3676 2578
21 4029 2838
22 4263 2931
23 3429 1908
24 3507 2367
25 3394 2643
26 3439 3201
27 2935 3240
28 3140 3550
29 2545 2357
30 2778 2826
31 2370 2975
以下為遺傳算法的主函數：
1 /*
2 文件名：CppGATSP 3 作者：Alex Xu 4 地址：Dalian Maritime University 5 描述：利用遺傳算法求解TSP問題(C++版) 6 創建時間：2018年12月10日11點27分 7 */
8
9 #include<iostream>
10 #include<vector>
11 #include<numeric> //accumulate
12 #include<chrono> //time
13 #include "HeuristicOperator.h"
14 using namespace std; 15 using namespace chrono; 16
17 //設置算法參數
18 # define POP_SIZE 2
19 # define MAX_GEN 4000
20
21 int main() { 22 //計時開始
23 auto start = system_clock::now(); 24
25 //生成距離矩陣
26 HeuristicOperator ga_dm; 27 vector<vector<double>> GA_DM; 28 GA_DM = ga_dm.getDM(ga_dm.getCoord()); 29
30 int n = int(GA_DM[0].size()); //城市規模 31
32 //初始化算法
33 vector<vector<int>> initPop(POP_SIZE, vector<int>(n)); //初始種群
34 vector<vector<int>> Pop(POP_SIZE, vector<int>(n)); //當前種群
35 vector<vector<int>> newPop(POP_SIZE, vector<int>(n)); //新種群
36 vector<double> popFit(POP_SIZE); //記錄種群適應度值
37 vector<int> bestIndival(n); //最優個體
38 vector<double> gs(MAX_GEN + 1); //記錄全局最優解
39 gs[0] = 1e9; 40 unsigned int seed = (unsigned)std::chrono::system_clock::now().time_since_epoch().count(); 41
42 //生成初始種群
43 HeuristicOperator s0; 44 for (int i = 0; i < POP_SIZE; i++) { 45 initPop[i] = s0.getInitS(n); 46 } 47 Pop = initPop; 48
49 //開始進化
50 for (int gen = 1; gen <= MAX_GEN; gen++) { 51
52 HeuristicOperator eval; //計算種群的適應度值（這里直接用路徑長度表示）
53 for (int i = 0; i < POP_SIZE; i++) { 54 popFit[i] = eval.Eval(Pop[i], GA_DM, n); 55 } 56
57 HeuristicOperator bestEI; //找出種群中個體的最優適應度值並記錄相應的個體編號
58 vector<double> bestEvalIndex(2); 59 bestEvalIndex = bestEI.bestS(popFit, POP_SIZE); 60 double bestEval = bestEvalIndex[0]; //最優適應度值
61 int bestIndex = int(bestEvalIndex[1]); //最優適應度值對應的個體編號 62
63 //最優解的更新
64 if (bestEval < gs[gen - 1]) { //比上一代優秀則更新
65 gs[gen] = bestEval; 66 bestIndival = Pop[bestIndex]; 67 } 68 else { //不比上一代優秀則不更新
69 gs[gen] = gs[gen - 1]; 70 } 71 if (gen % 100 == 0) { 72 cout << "第" << gen << "次迭代時全局最優評價值為" << gs[gen] << endl; 73 } 74
75 //擾動操作（產生新種群）
76 for (int p = 0; p < POP_SIZE; p++) { 77 HeuristicOperator shk; 78 vector<int> randPosition = shk.RandPosition(n); 79 vector<int> tmpS(n); 80 double randShk = rand() / double(RAND_MAX); 81 if (randShk < 0.33) { 82 tmpS = shk.Swap(Pop[p], randPosition); //交換操作
83 } 84 else if (randShk >= 0.67) { 85 tmpS = shk.Flip(Pop[p], randPosition); //翻轉操作
86 } 87 else { 88 tmpS = shk.Insert(Pop[p], randPosition); //插入操作
89 } 90
91 HeuristicOperator evl; 92 if (evl.Eval(tmpS, GA_DM, n) > evl.Eval(Pop[p], GA_DM, n)) { 93 newPop[p] = Pop[p]; 94 } 95 else { 96 newPop[p] = tmpS; 97 } 98 } 99 Pop = newPop; 100
101 //選擇操作（輪盤賭）
102 vector<double> Cusum(POP_SIZE + 1, 0); //適用於輪盤賭的累加器Cusum（補充了cus[0]=0;
103 for (int i = 0; i < POP_SIZE; i++) { 104 Cusum[i + 1] = Cusum[i] + popFit[i]; 105 } 106
107 double Sum = accumulate(popFit.begin(), popFit.end(), 0.0); //計算各個體被選擇的概率（歸一化）
108 vector<double> cusFit(POP_SIZE + 1); //放置種群中個個體被選擇的概率值
109 for (int i = 0; i < POP_SIZE + 1; i++) { 110 cusFit[i] = Cusum[i] / Sum; 111 } 112
113 for (int p = 0; p < POP_SIZE; p++) { //輪盤賭產生新種群
114 double r = rand() / double(RAND_MAX); 115 for (int q = 0; q < POP_SIZE; q++) { 116 if (r > cusFit[q] && r <= cusFit[q + 1]) { 117 newPop[p] = Pop[q]; 118 } 119 } 120 } 121 Pop = newPop; 122 } 123
124 //計時結束
125 auto end = system_clock::now(); 126 auto duration = duration_cast<microseconds>(end - start); 127 cout << "花費了"
128 << double(duration.count()) * microseconds::period::num / microseconds::period::den 129 << "秒" << endl; 130
131 //輸出結果
132 double gs0 = 15377.711; 133 cout << "最優解為" << gs[MAX_GEN] << endl; 134 double e = (gs[MAX_GEN] - gs0) / gs0; 135 cout << "誤差為" << e * 100.0 << '%' << endl; 136 cout << "最優路徑為" << endl; 137 for (int i = 0; i < n; i++) { 138 cout << bestIndival[i] + 1 << '\t'; 139 }
140
141 while (1)
142 {}
143 }
以上即為C++語言所編寫的遺傳算法求解TSP示例，運行環境為：Windows10 64位操作系統；CPU：i7-8750H; 內存：8G；Microsoft Visual Studio Community 2017 。求解結果如下：
與已知最優解的誤差為1.48%，所用時間約為3.6s. 還可以接受。但值得注意的是：本文實驗參數最大迭代次數4000代，而種群規模僅為2，這與一般的遺傳算法思想上是沒問題的，只是實際參數可能取得不太大眾化。當然，對於算法參數這些細節都是可以調節的，不必太過於糾結。
啊哈，這次的利用C++編程遺傳算法求解TSP就這些了~~~