用遺傳算法解決旅行商問題(附源代碼)
本文程序所有源代碼已在《用遺傳算法解決旅行商問題開源(全部源代碼)》中開源。
最近心血來潮，重新拾起大學畢業設計時研究過的遺傳算法。去年做畢業設計時還覺得遺傳算法是一種多么神秘的算法，但是今天看來，遺傳算法也就和冒泡排序算法差不多，都是通用的算法，只不過遺傳算法實現起來稍微復雜一點而已。
我曾經被遺傳算法的名字所疑惑，還以為遺傳算法會改變程序的形態，使得程序就好像生物一樣進化，過了幾天去看程序已經變得連編寫程序的人都認不出來了，汗！大二時的幼稚想法。
遺傳算法其實是一種求函數極值的隨機搜索算法，但它又不是毫無規則地隨機搜索，而是基於一種假設：假設函數值的分布是有一定的連續性的，換句話說函數的極值出現在一個較優值附近的概率要大於出現在一個較差值附近的概率。基於這個假設，遺傳算法總是以較大概率保留較優值所代表的搜索方向，而以較低概率保留較差值所代表的搜索方向。這並不是說不去搜索較差值的附近區域，只是搜索的概率較低而已。這個思想與模擬退火算法相似，對於能量較高的系統狀態，程序仍然以一定的概率接受，只不過這個概率小於1。
遺傳算法的局部搜索能力較強，但是很容易陷入局部極值，畢業設計的時候曾經認為只要增加變異概率就可以跳出局部極值，還美其名曰自適應，現在想想這種想法是錯誤的：雖然增加變異概率可以搜索到遠離當前極值的點，但是新點的值往往不能和當前保留下來的較優值相提並論，因為這些較優值都是經過千百代的進化而存留下來的，於是遠離當前極值的點往往在兩到三代以內就被淘汰掉了。增加變異概率實際上是把遺傳算法退化成了一種純粹的隨機搜索，所謂的自適應也無從談起！
那么如何解決遺傳算法容易陷入局部極值的問題呢？讓我們來看看大自然提供的方案。六千五百萬年以前，恐龍和靈長類動物並存，恐龍在地球上占絕對統治地位，如果恐龍沒有滅絕靈長類動物是絕沒有可能統治地球的。正是恐龍的滅絕才使靈長類動物有了充分進化的余地，事實上地球至少經歷了5次物種大滅絕，每次物種滅絕都給更加高級的生物提供了充分進化的余地。所以要跳出局部極值就必須殺死當前所有的優秀個體，從而讓遠離當前極值的點有充分的進化余地。這就是災變的思想。
下一個問題是什么時候進行災變，換句話說什么時候局部搜索已經充分了呢？我用了一個災變倒計數的概念：從500開始遞減，每一代遞減一次，如果出現了新的最優值，就從新開始計數，如果出現新最優值的時候倒計數遞減次數的2.5倍已經超過500則從新的初始值開始倒數。例：初始倒數500，如果倒數到200時出現新最優值，則從(500 - 200) * 2.5 = 750開始從新倒數，下一次如果倒數到100時出現新最優值，則從(750 - 100) * 2.5 = 1625開始倒計數(這里的2.5是一個經驗值，可以在全局參數設置里面調整)。也就是說倒計數的長度取決於進化的速度，進化速度越慢倒計數長度越長。如果倒計數完畢還沒有新的最優值，就認為局部搜索已經充分，就發生災變。
基於上訴思想我寫了一個程序來計算旅行商問題。我現在終於體會到旅行商問題為什么會這么有名，有很多算法都可以解決旅行商問題，問題描述簡單，評價函數也不復雜，問題的解可以直觀地顯示出來，具有各種如局部極值多等典型的性質，這些都成為算法練兵的好處，可以清晰地比較各個算法的優劣，發現算法的缺陷。可以說旅行商問題就是一個練兵場，一個學校，為算法提供了成長的場所。為算法能夠應用到其他復雜領域打好基礎。
程序輸入是一個文本文件，里面記錄了所有城市的坐標，以及最優個體的序列。以一張只有10個城市的地圖為例，文本中可能記錄了以下內容：
0.604600, 0.592500, 8
0.610500, 0.261400, 3
0.572800, 0.494300, 7
0.153200, 0.983900, 2
0.955700, 0.772000, 0
0.914400, 0.276500, 4
0.998500, 0.484800, 6
0.449800, 0.605300, 5
0.308500, 0.109000, 1
0.364700, 0.060100, 9
表示第一個城市的坐標為0.604600, 0.592500(程序客戶區的寬和高為單位1，所有城市的坐標值均在[0.0,0.0] ~ [1.0,1.0]之內)，第二個城市坐標為0.610500, 0.261400...依次類推。
后面所跟的整數為最優個體的序列，上述數據表示旅行商應該從第8號城市(0.308500, 0.109000)出發，經過3，7，2，0，4，6，5，1，9號城市，最后又回到第8號城市。
程序的最終目標是求取一個序列，使得旅行商按照這個序列旅行時行程最短。
程序的變異方法是自繁殖變異，有兩種：1、隨機取兩點，逆序這兩點間的序列。2、隨機把一個城市轉移到另一個序列位置。
對於一個500個城市的地圖，大概在5萬代左右發生第一次災變，用時約6~8分鍾，災變前夕的災變倒計數初始值已經從800達到2000~20000。可以看到從一次災變結束到下一次災變開始，最優值的變化趨勢近似呈一條拖拽線，越接近局部極值進化速度越慢，這也說明災變倒計數的策略是正確的。
下面是一次試驗的數據統計：程序運行兩個小時，進化到一百萬代，發生了16次災變，最優個體產生於第606722代，屬於第11個進化周期，總行程長度為17.164006，第一次災變發生在第49773代，第一次災變前最優個體產生於第45523代，總行程長度為18.029128。
圖1 最佳路線圖
圖2 第一次災變前的最佳路線
分析最優分趨勢圖
圖3 第一次災變前的最優分趨勢圖
圖4 最后一次災變前的最優分趨勢圖
在每個進化周期內最優分圖形基本呈拖拽線形狀，可以看到多數進化周期已經沒有進化速度，說明局部搜索已經充分，少數進化周期在發生災變時還有明顯進化速度，這是因為這些周期恰好進入一個比較長的停滯期時被程序認為局部搜索充分了，停滯期的出現根隨機數有關，個人認為應該可以通過調節災變初始值和災變倍增值解決。
分析平均分趨勢圖
圖5 第一次災變前的平均分趨勢圖
圖6 最后一次災變前的平均分趨勢圖
可以看到每次分生災變后群體平均分會達到一個較大的值，然后迅速下降，再慢慢上升。這說明旅行商問題的局部極值非常多，極值附近解的分數要遠遠低於整個解空間的平均分，這主要是因為一個較優解的進行一次變異后生成的子女絕大部分都是畸形的分數很低的個體，由於遺傳算法並不放棄這些進化方向，從而影響了群體的平均分。災變時對整個解空間進行隨機搜索，這時的群體平均分可以作為整個解空間平均分的體現，進化一定時間以后，群體迅速陷入到一個特定的局部極值附近，這個時候較優解還沒有進化出來，群體中充斥着畸形個體，只有少量比較優秀的個體，所以平均分也隨之迅速下降，隨后由於優秀個體存活率比較高，群體漸漸被優秀個體統治，群體平均分也開始上升。仔細分析每一個進化周期的平均分趨勢圖可以發現，在進化的后期群體平均分有一個穩固上升的階段（這應該是最優個體慢慢排擠其他個體的結果），在此之前都會有一個標志性的少量下挫曲線（如圖），還不知道產生這個曲線的原因。
理論上說，保留更多的種子，可以保留更多的搜索方向，搜索的效果應該越好，但是我做了一下對比試驗，發現保留1個種子的搜索速度更快，搜索到的極值更優秀，但是這與遺傳算法的精髓相違背，困惑中。
程序下載地址：
下面是程序的主要核心代碼：
// 變異函數
inline void Variant(GENE & gsource, GENE & gdest, int * ptemp, int size, int varate)
{
static int i;
static int j, k, l, n, m;
static int tmp;
memcpy(gdest.index, gsource.index, sizeof(int) * size);
for(i = 0; i < varate; i++)
{
switch(rand() % 2)
{
case 0:
// 逆序變異
{
j = rand() % size;
k = rand() % size;
if(j == k)
{
k = (k + 1) % size;
}
if(j > k)
{
k += size;
for(l = j; l < j + k - l; l++)
{
n = (j + k - l) % size;
m = l % size;
tmp = gdest.index[m];
gdest.index[m] = gdest.index[n];
gdest.index[n] = tmp;
}
}
else
{
for(l = j; l < j + k - l; l++)
{
tmp = gdest.index[l];
gdest.index[l] = gdest.index[j + k - l];
gdest.index[j + k - l] = tmp;
}
}
}
break;
case 1:
// 轉移變異
{
j = rand() % size;
k = 1;//rand() % ((size >> 4) + 1);
if(k == 0)
{
k = 1;
}
l = rand() % (size - k) + 1;
n = (j + k) % size + l > size ? size - (j + k) % size : l;
memcpy(ptemp, gdest.index + (j + k) % size, n * sizeof(int));
if(n < l)
{
memcpy(ptemp + n, gdest.index, (l - n) * sizeof(int));
}
n = j + k > size ? size - j : k;
memcpy(ptemp + l, gdest.index + j, n * sizeof(int));
if(n < k)
{
memcpy(ptemp + l + n, gdest.index, (k - n) * sizeof(int));
}
m = k + l;
n = j + m > size ? size - j :m;
memcpy(gdest.index + j, ptemp, n * sizeof(int));
if(n < m)
{
memcpy(gdest.index, ptemp + n, (m - n) * sizeof(int));
}
}
break;
default:
break;
}
}
}
// 輔助線程
UINT CCityMap::ThreadProc( LPVOID pParam )
{
CCityMap * pClass = (CCityMap *)pParam;
if(pClass->m_iCityNum <= 1)
{
return 0;
}
srand((UINT)time(NULL));
pClass->m_bCompute = true;
pClass->ReadSetting();
static int i, j;
static GENE tgene;
static int bullet;
static int totalBullet;
static int maxCountdown = CM_JUMP_COUNTDOWN_INIT;
// 分配空間
int num = pClass->m_iCityNum;
int * ptmp = new int[num];
int commSize = CM_SEED_NUM * (1 + CM_CHILDREN_NUM);
GENE * comm = new GENE[commSize];
for(i = 0; i < commSize; i++)
{
comm[i].index = new int [num];
comm[i].mark = 0.0;
}
// 生成初始群落
CTime tstart = CTime::GetCurrentTime();
CTime tnow;
pClass->m_i64GenNum = 1;
pClass->m_iJumpCount = 0;
memcpy(comm[0].index, pClass->m_piBestIndex, sizeof(int) * num);
pClass->Mark(comm[0]);
pClass->QuadrangleOptimise(comm[0]);
pClass->m_dBestMark = comm[0].mark;
pClass->m_i64BestGen = pClass->m_i64GenNum;
double maxMark = comm[0].mark;
int maxIndex = 0;
double totalMark = maxMark;
pClass->m_iJumpCountdown = CM_JUMP_COUNTDOWN_INIT;
for(i = 1; i < CM_SEED_NUM; i++)
{
Variant(comm[0], comm[i], ptmp, num, pClass->m_iJumpCountdown);
pClass->Mark(comm[i]);
totalMark += comm[i].mark;
if(maxMark < comm[i].mark)
{
maxMark = comm[i].mark;
maxIndex = i;
}
}
// 移動最優基因
if(maxIndex != 0)
{
tgene.index = comm[0].index;
tgene.mark = comm[0].mark;
comm[0].index = comm[maxIndex].index;
comm[0].mark = comm[maxIndex].mark;
comm[maxIndex].index = tgene.index;
comm[maxIndex].mark = tgene.mark;
maxIndex = 0;
pClass->QuadrangleOptimise(comm[0]);
maxMark = comm[0].mark;
memcpy(pClass->m_piBestIndex, comm[0].index, sizeof(int) * num);
pClass->m_dBestMark = maxMark;
pClass->m_i64BestGen = pClass->m_i64GenNum;
}
int indNum = CM_SEED_NUM;
pClass->m_dAVGMark = (totalMark) / indNum;
pClass->m_pMaxTrendLine->Clear();
pClass->m_pAVGTrendLine->Clear();
pClass->m_pMaxTrendLine->AddValue(maxMark);
pClass->m_pAVGTrendLine->AddValue(pClass->m_dAVGMark);
// 開始進化
while(!pClass->m_bKillMsg)
{
totalMark = 0.0;
// 變異
for(i = 0; i < CM_SEED_NUM; i++)
{
totalMark += comm[i].mark;
for(j = 0; j < CM_CHILDREN_NUM; j++)
{
Variant(comm[i], comm[indNum], ptmp, num, 1);
pClass->Mark(comm[indNum]);
totalMark += comm[indNum].mark;
if(maxMark < comm[indNum].mark)
{
maxMark = comm[indNum].mark;
maxIndex = indNum;
}
indNum++;
}
}
pClass->m_dAVGMark = (totalMark) / indNum;
pClass->m_pAVGTrendLine->AddValue(pClass->m_dAVGMark);
pClass->m_pMaxTrendLine->AddValue(maxMark);
// 移動最優基因
if(maxIndex != 0)
{
tgene.index = comm[0].index;
tgene.mark = comm[0].mark;
comm[0].index = comm[maxIndex].index;
comm[0].mark = comm[maxIndex].mark;
comm[maxIndex].index = tgene.index;
comm[maxIndex].mark = tgene.mark;
maxIndex = 0;
if(maxMark > pClass->m_dBestMark)
{
memcpy(pClass->m_piBestIndex, comm[0].index, sizeof(int) * num);
pClass->m_dBestMark = maxMark;
pClass->m_i64BestGen = pClass->m_i64GenNum;
}
int forcastCountdown = int((maxCountdown - pClass->m_iJumpCountdown) * CM_JUMP_COUNTDOWN_INC);
if(forcastCountdown > maxCountdown)
{
maxCountdown = forcastCountdown;
}
pClass->m_iJumpCountdown = maxCountdown;
}
else
{
pClass->m_iJumpCountdown--;
if(pClass->m_iJumpCountdown <= 0)
{
pClass->QuadrangleOptimise(comm[0]);
if(maxMark < comm[0].mark)
{
pClass->m_iJumpCountdown = maxCountdown;
maxMark = comm[0].mark;
if(maxMark > pClass->m_dBestMark)
{
memcpy(pClass->m_piBestIndex, comm[0].index, sizeof(int) * num);
pClass->m_dBestMark = maxMark;
pClass->m_i64BestGen = pClass->m_i64GenNum;
}
}
else
{
if(CM_IMG_LOG)
{
// 保存當前屏幕圖像為文件,作為日志
pClass->m_iJumpCountdown = maxCountdown;
static CString fileName;
fileName.Format("%03d.bmp", pClass->m_iJumpCount);
pClass->SaveAsImage(fileName);
}
srand((UINT)time(NULL));
maxCountdown = CM_JUMP_COUNTDOWN_INIT;
pClass->m_iJumpCountdown = maxCountdown;
// 已經陷入局部最優，災變
pClass->m_iJumpCount++;
Variant(comm[0], comm[0], ptmp, num, 20);
pClass->Mark(comm[0]);
maxMark = comm[0].mark;
for(i = 1; i < CM_SEED_NUM; i++)
{
Variant(comm[0], comm[i], ptmp, num, 20);
pClass->Mark(comm[i]);
totalMark += comm[i].mark;
if(maxMark < comm[i].mark)
{
maxMark = comm[i].mark;
maxIndex = i;
}
}
// 移動最優基因
if(maxIndex != 0)
{
tgene.index = comm[0].index;
tgene.mark = comm[0].mark;
comm[0].index = comm[maxIndex].index;
comm[0].mark = comm[maxIndex].mark;
comm[maxIndex].index = tgene.index;
comm[maxIndex].mark = tgene.mark;
maxIndex = 0;
}
indNum = CM_SEED_NUM;
}
}
}
// 輪盤賭
totalMark -= comm[0].mark;
totalBullet = 0;
for(i = 1; i < indNum; i++)
{
comm[i].killRate = int(10000.0 * comm[i].mark / totalMark);
totalBullet += comm[i].killRate;
}
while(indNum > CM_SEED_NUM)
{
bullet = rand() % totalBullet;
for(i = 1; i < indNum; i++)
{
if(bullet <= comm[i].killRate)
{
// 命中
totalBullet -= comm[i].killRate;
tgene.index = comm[indNum - 1].index;
tgene.mark = comm[indNum - 1].mark;
tgene.killRate = comm[indNum - 1].killRate;
comm[indNum - 1].index = comm[i].index;
comm[indNum - 1].mark = comm[i].mark;
comm[indNum - 1].killRate = comm[i].killRate;
comm[i].index = tgene.index;
comm[i].mark = tgene.mark;
comm[i].killRate = tgene.killRate;
indNum--;
break;
}
else
{
bullet -= comm[i].killRate;
}
}
}
pClass->m_i64GenNum++;
tnow = CTime::GetCurrentTime();
pClass->m_tsTimeUsed = tnow - tstart;
}
// 釋放空間
for(i = 0; i < commSize; i++)
{
delete [] comm[i].index;
comm[i].index = NULL;
}
delete [] comm;
comm = NULL;
tgene.index = NULL;
pClass->m_bCompute = false;
return 0; // thread completed successfully
}