模擬退火演算法與C語言實現(TSP問題)
1簡介:
模擬退火來自冶金學的專有名詞退火。退火是將材料加熱後再經特定速率冷卻,目的是增大晶粒的體積,並且減少晶格中的缺陷。材料中的原子原來會停留在使內能有區域性最小值的位置,加熱使能量變大,原子會離開原來位置,而隨機在其他位置中移動。退火冷卻時速度較慢,使得原子有較多可能可以找到內能比原先更低的位置。
模擬退火的原理也和金屬退火的原理近似:我們將熱力學的理論套用到統計學上,將搜尋空間內每一點想象成空氣內的分子;分子的能量,就是它本身的動能;而搜尋空間內的每一點,也像空氣分子一樣帶有“能量”,以表示該點對命題的合適程度。演算法會計算原有的解的適合度與新的解的適合度的差值,通過Metropolis準則計算出概率來決定是否接受新的解。在Metropolis準則之下,在開始的階段更容易接受較差的解,而隨著時間的推移,較差的解越來越難以被接受,並最終穩定(即收斂)。這種接受較差的解的特性,使得模擬退火演算法有能力跳出區域性最優解,找到全域性最優解。
模擬退火是一種通用概率演算法,常用來在一定時間內尋找在一個很大搜尋空間中的近似最優解。
2演算法要點:
從上述簡介中不難總結出演算法設計的幾個要點,首先要設計檢測命題適合程度的函式,從而計算適合度的大小進行比較;另外在進行計算概率的時候根據的是Metropolis準則,即在溫度T時趨於平衡的概率為,其中E為溫度T時的內能,ΔE為其改變數,k為Boltzmann常數。其中的ΔE在演算法中化為兩個解適合度的差值,T為控制引數。在T足夠大的時候,適合度差值較大的情況也能有較大的機率接受新較差的解,而隨著T不斷變小,適合度差值也只有在較小的範圍內才能接受新的較差的解,並最終穩定,輸出最優解。
3
初始化:
定義初始溫度引數T和降溫係數q,隨機生成一組當前問題的解,根據適合度函式計算適度值f1。
迭代過程:
迭代過程是模擬退火演算法的核心步驟,分為新解的產生和接受新解兩部分:
1)由一個產生函式從當前解產生一個位於解空間的新解;為便於後續的計算和接受,減少演算法耗時,通常選擇由當前新解經過簡單地變換即可產生新解的方法,如對構成新解的全部或部分元素進行置換、互換等,注意到產生新解的變換方法決定了當前新解的鄰域結構,因而對冷卻進度表的選取有一定的影響。
2)計算與新解所對應的目標函式差。因為目標函式差僅由變換部分產生,所以目標函式差的計算最好按增量計算。事實表明,對大多數應用而言,這是計算目標函式差的最快方法。
3)判斷新解是否被接受,判斷的依據是一個接受準則,最常用的接受準則是Metropolis準則:若Δt′<0則接受S′作為新的當前解S,否則以概率exp(-Δt′/T)接受S′作為新的當前解S。
4)當新解被確定接受時,用新解代替當前解,這隻需將當前解中對應於產生新解時的變換部分予以實現,同時修正目標函式值即可。此時,當前解實現了一次迭代。可在此基礎上開始下一輪試驗。而當新解被判定為捨棄時,則在原當前解的基礎上繼續下一輪試驗。
模擬退火演算法與初始值無關,演算法求得的解與初始解狀態S(是演算法迭代的起點)無關;模擬退火演算法具有漸近收斂性,已在理論上被證明是一種以概率1收斂於全域性最優解的全域性優化演算法;模擬退火演算法具有並行性。
停止準則:
溫度T降至某最低值時,完成給定數量迭代中無法接受新解,停止迭代,接受當前尋找的最優解為最終解。
退火方案:
在某個溫度狀態T下,當一定數量的迭代操作完成後,降低溫度T,在新的溫度狀態下執行下一個批次的迭代操作。
程式碼如下:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#define T_start 5000.0 //初始溫度
#define T_end (1e-8) //結束溫度
#define q 0.98 //退火係數
#define L 1000 //每個溫度最大迭代次數
#define N 34 //城市個數
int city_result[N]; //城市列表的解空間
double city[N][2] = {{9932, 4439}, {10109, 4351}, {11552, 3472}, {10302, 3290}, {8776, 3333}, {7040, 4867}, {9252, 4278}, {9395, 4539}, {11101, 2540}, {9825, 5087}, {10047, 4879}, {10227, 4648}, {100027, 4229}, {9878, 4211}, {9087, 4065}, {10438, 4075}, {10382, 3865}, {11196, 3563}, {11075, 3543}, {11544, 3365}, {11915, 2900}, {11305, 3189}, {11073, 3137}, {10950, 3394}, {11576, 2575}, {12239, 2785}, {11529, 2226}, {9328, 4006}, {10012, 3811}, {9952, 3410}, {10612, 2954}, {10349, 2784}, {11747, 2469}, {11673, 2461}}; //中國34個城市實際距離座標
//函式宣告
//double city[N][2] = {{0,0},{3,0},{3,4}};//測試資料
double distance(double *city1, double *city2); //計算兩城市間距離
double path(int city_result[N]); //計算總路徑
void init(); //初始化解空間
void creat(); //生成新的解空間
int main()
{
time_t start, end; //記錄程式開始結束時間
double time_sum; //記錄程式執行時間
start = clock(); //程式開始時間
int i, count; //降溫計數器
int city_copyresult[N]; //拷貝解空間
double path1, path2; //原有解空間,新解空間的總路徑
double dE; //原有解空間與新解空間的差值
double r; //隨機產生0~1的值,是否接受新的解
double T; //當前溫度
srand(time(NULL));
init(); //初始化解空間
T = T_start;//初始溫度賦值
while (T > T_end) //當前溫度大於結束溫度
{
for (i = 0; i < L; i++)
{
memcpy(city_copyresult, city_result, N * sizeof(int));
creat(); //產生新的解空間
path1 = path(city_copyresult);
path2 = path(city_result);
dE = path2 - path1;
if (dE > 0) //Metropolis準則
{
r = rand() / (RAND_MAX);
if (exp(-dE / T) <= r) //高溫狀態下,可以接受能量差值較大的新狀態;低溫狀態下,則只能接受能量差值較小的新狀態
memcpy(city_result, city_copyresult, N * sizeof(int)); //保留原來的解
}
}
T *= q;
count++;
}
end = clock(); //程式結束時間
time_sum = (double)(end - start ) / (CLOCKS_PER_SEC); //換算成秒
printf("共降溫:%d次\n", count);
printf("經過模擬退火演算法得出最優路徑長度為:%f\n", path(city_result));
for (i = 0; i < N; i++)
{
printf("%d->", city_result[i]);
}
printf("%d\n", city_result[0]);
printf("程式共耗時%f秒.\n",time_sum);
system("pause");
return 0;
}
double distance(double *city1, double *city2) //計算兩個城市之間的距離
{
double x1, x2, y1, y2, dis;
x1 = *city1; //第一個城市的x座標
x2 = *city2; //第二個城市的x座標
y1 = *(city1 + 1); //第一個城市的y座標
y2 = *(city2 + 1); //第二個城市的y座標
dis = sqrt(pow((x2 - x1), 2) + pow((y2 - y1), 2));
return dis;
}
double path(int city_result[N]) //計算總路徑
{
int i, city1_num, city2_num; //解空間中的兩個城市序號
double sum = 0; //路徑總長度
for (i = 0; i < N - 1; i++) //解空間中首位到末位的總路徑
{
city1_num = city_result[i];
city2_num = city_result[i + 1];
sum += distance(city[city1_num], city[city2_num]);
}
sum += distance(city[0], city[N - 1]); //加上解空間中末位到首位的路徑
return sum;
}
void init() //初始化解空間
{
int i;
for (i = 0; i < N; i++)
city_result[i] = i; //順序生成解空間
}
void creat() //生成新的解空間
{
int point1, point2, temp;
point1 = rand() % N;
point2 = rand() % N;
temp = city_result[point1];
city_result[point1] = city_result[point2];
city_result[point2] = temp;
}
本文參考:
https://zh.wikiphttps://zh.wikipedia.org/wiki/%E6%A8%A1%E6%8B%9F%E9%80%80%E7%81%ABedia.org/wiki/%E6%A8%A1%E6%8B%9F%E9%80%80%E7%81%AB
https://www.cnblogs.com/lyrichu/p/6688459.html