C++版遺傳演算法求解TSP
隔半年，再次有時間整理關於組合優化問題——旅行商問題（Traveling Salesman Problem, TSP），這次採用的是經典遺傳演算法（Genetic Algorithm， GA）進行求解，利用C++語言進行程式設計實現。
各種啟發式演算法的整體框架大致都由以下幾個操作組成：（1）初始解的產生；（2）解的評價（評價函式）；（3）擾動運算元；此外，還可以加上程式原始資料的匯入等操作。這些操作是多數啟發式演算法所通用的運算元，基為此，此次在採用C++進行實現的時候，採用一個通用的 HeuristicOperator.h 標頭檔案以及對應的 HeuristicOperator.cpp 類檔案對這些操作進行集中放置，造好輪子，方便以後取用。
表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標頭檔案：
#pragma once
#include<iostream>
#include<vector>
#include <algorithm> // std::shuffle
#include <random> // std::default_random_engine
#include<chrono>
using namespace std;
class HeuristicOperator {
public:
vector<vector<double>> getCoord(void); //函式1：獲取座標函式
vector<vector<double>> getDM(vector<vector<double>> Coord); //函式2：獲取距離矩陣函式
vector<int> getInitS(int n); //函式3：獲取初始解函式
double Eval(vector<int> S, vector<vector<double>> DM, int n); //函式4：評價函式
vector<double> bestS(vector<double> Eval, int Length); //函式5：搜尋範圍內最優評價值及其相應的位置函式
vector<int> RandPosition(int n); //函式6：產生Sharking操作位置函式
vector<int> Swap(vector<int> S, vector<int> RP); //函式7：交換運算元
vector<int> Flip(vector<int> S, vector<int> RP); //函式8：翻轉運算元
vector<int> Insert(vector<int> S, vector<int> RP); //函式9：插入運算元
};
對應的HeuristicOperator.cpp類檔案比較容易實現，在此不再贅述。本文所用算例為31城市的TSP問題，與Java版遺傳演算法求解TSP求解算例一致，具體資料如下：
1304 2312
3639 1315
4177 2244
3712 1399
3488 1535
3326 1556
3238 1229
4196 1004
4312 790
4386 570
3007 1970
2562 1756
2788 1491
2381 1676
1332 695
3715 1678
3918 2179
4061 2370
3780 2212
3676 2578
4029 2838
4263 2931
3429 1908
3507 2367
3394 2643
3439 3201
2935 3240
3140 3550
2545 2357
2778 2826
2370 2975
以下為遺傳演算法的主函式：
/*
檔名：CppGATSP
作者：Alex Xu
地址：Dalian Maritime University
描述：利用遺傳演算法求解TSP問題
建立時間：2018年12月10日11點27分
*/
#include<iostream>
#include<vector>
#include<numeric> //accumulate
#include<chrono> //time
#include "HeuristicOperator.h"
using namespace std;
using namespace chrono;
//設定演算法引數
# define POP_SIZE 2
# define MAX_GEN 4000
int main() {
//計時開始
auto start = system_clock::now();
//生成距離矩陣
HeuristicOperator ga_dm;
vector<vector<double>> GA_DM;
GA_DM = ga_dm.getDM(ga_dm.getCoord());
int n = int(GA_DM[0].size()); //城市規模
//初始化演算法
vector<vector<int>> initPop(POP_SIZE, vector<int>(n)); //初始種群
vector<vector<int>> Pop(POP_SIZE, vector<int>(n)); //當前種群
vector<vector<int>> newPop(POP_SIZE, vector<int>(n)); //新種群
vector<double> popFit(POP_SIZE); //記錄種群適應度值
vector<int> bestIndival(n); //最優個體
vector<double> gs(MAX_GEN + 1); //記錄全域性最優解
gs[0] = 1e9;
unsigned int seed = (unsigned)std::chrono::system_clock::now().time_since_epoch().count();
//生成初始種群
HeuristicOperator s0;
for (int i = 0; i < POP_SIZE; i++) {
initPop[i] = s0.getInitS(n);
}
Pop = initPop;
//開始進化
for (int gen = 1; gen <= MAX_GEN; gen++) {
HeuristicOperator eval; //計算種群的適應度值（這裡直接用路徑長度表示）
for (int i = 0; i < POP_SIZE; i++) {
popFit[i] = eval.Eval(Pop[i], GA_DM, n);
}
HeuristicOperator bestEI; //找出種群中個體的最優適應度值並記錄相應的個體編號
vector<double> bestEvalIndex(2);
bestEvalIndex = bestEI.bestS(popFit, POP_SIZE);
double bestEval = bestEvalIndex[0]; //最優適應度值
int bestIndex = int(bestEvalIndex[1]); //最優適應度值對應的個體編號
//最優解的更新
if (bestEval < gs[gen - 1]) { //比上一代優秀則更新
gs[gen] = bestEval;
bestIndival = Pop[bestIndex];
}
else { //不比上一代優秀則不更新
gs[gen] = gs[gen - 1];
}
if (gen % 100 == 0) {
cout << "第" << gen << "次迭代時全域性最優評價值為" << gs[gen] << endl;
}
//擾動操作（產生新種群）
for (int p = 0; p < POP_SIZE; p++) {
HeuristicOperator shk;
vector<int> randPosition = shk.RandPosition(n);
vector<int> tmpS(n);
double randShk = rand() / double(RAND_MAX);
if (randShk < 0.33) {
tmpS = shk.Swap(Pop[p], randPosition); //交換操作
}
else if (randShk >= 0.67) {
tmpS = shk.Flip(Pop[p], randPosition); //翻轉操作
}
else {
tmpS = shk.Insert(Pop[p], randPosition); //插入操作
}
HeuristicOperator evl;
if (evl.Eval(tmpS, GA_DM, n) > evl.Eval(Pop[p], GA_DM, n)) {
newPop[p] = Pop[p];
}
else {
newPop[p] = tmpS;
}
}
Pop = newPop;
//選擇操作（輪盤賭）
vector<double> Cusum(POP_SIZE + 1, 0); //適用於輪盤賭的累加器Cusum（補充了cus[0]=0;
for (int i = 0; i < POP_SIZE; i++) {
Cusum[i + 1] = Cusum[i] + popFit[i];
}
double Sum = accumulate(popFit.begin(), popFit.end(), 0.0); //計算各個體被選擇的概率（歸一化）
vector<double> cusFit(POP_SIZE + 1); //放置種群中個個體被選擇的概率值
for (int i = 0; i < POP_SIZE + 1; i++) {
cusFit[i] = Cusum[i] / Sum;
}
for (int p = 0; p < POP_SIZE; p++) { //輪盤賭產生新種群
double r = rand() / double(RAND_MAX);
for (int q = 0; q < POP_SIZE; q++) {
if (r > cusFit[q] && r <= cusFit[q + 1]) {
newPop[p] = Pop[q];
}
}
}
Pop = newPop;
}
//計時結束
auto end = system_clock::now();
auto duration = duration_cast<microseconds>(end - start);
cout << "花費了"
<< double(duration.count()) * microseconds::period::num / microseconds::period::den
<< "秒" << endl;
//輸出結果
double gs0 = 15377.711;
cout << "最優解為" << gs[MAX_GEN] << endl;
double e = (gs[MAX_GEN] - gs0) / gs0;
cout << "誤差為" << e * 100.0 << '%' << endl;
cout << "最優路徑為" << endl;
for (int i = 0; i < n; i++) {
cout << bestIndival[i] + 1 << '\t';
}
while (1)
{}
}
以上即為C++語言所編寫的遺傳演算法求解TSP示例，執行環境為：Windows10 64位作業系統；CPU：i7-8750H; 記憶體：8G；Microsoft Visual Studio Community 2017 。求解結果如下：
與已知最優解的誤差為1.48%，所用時間約為3.6s. 還可以接受。但值得注意的是：本文實驗引數最大迭代次數4000代，而種群規模僅為2，這與一般的遺傳演算法思想上是沒問題的，只是實際引數可能不太符合。當然，對於演算法引數這些細節都是可以調節的，不必太過於糾結。
啊哈，這次的利用C++程式設計遺傳演算法求解TSP就這些了~~~