1. 程式人生 > >模擬退火演算法與C語言實現(TSP問題)

模擬退火演算法與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