天天看點

數學模組化暑期集訓23:模拟退火算法模拟退火算法類屬爬山法模拟退火算法流程 C t {C_t} Ct​的設定引入 C t {C_t} Ct​後的模拟退火算法流程例題1:求給定函數的最大值或者最小值例題2:旅行商問題總結

模拟退火算法類屬

和粒子群算法一樣,模拟退火算法也屬于啟發式算法的一種。

啟發式算法,可參照下面的定義。

啟發式算法:在搜尋最優解的過程中利用到了原來搜尋過程中得到的資訊,且這個資訊會改進我們的搜尋過程。

爬山法

模拟退火算法,可以算一種優化過的爬山法。

數學模組化暑期集訓23:模拟退火算法模拟退火算法類屬爬山法模拟退火算法流程 C t {C_t} Ct​的設定引入 C t {C_t} Ct​後的模拟退火算法流程例題1:求給定函數的最大值或者最小值例題2:旅行商問題總結

爬山法比較好了解,首先在函數圖上随機選取一個點,之後再其左邊或右邊各選一點,若比該點大,則以大的點繼續選擇,整個過程類似于爬山。

問題在于,當爬到小山峰的時候,無法繼續爬,這就導緻陷入局部最優解。

模拟退火算法流程

模拟退火在爬山法的基礎上,結合蒙特卡洛的思想,整個流程如下:

數學模組化暑期集訓23:模拟退火算法模拟退火算法類屬爬山法模拟退火算法流程 C t {C_t} Ct​的設定引入 C t {C_t} Ct​後的模拟退火算法流程例題1:求給定函數的最大值或者最小值例題2:旅行商問題總結

C t {C_t} Ct​的設定

在上面的算法流程中, C t {C_t} Ct​并沒有具體設定, C t {C_t} Ct​控制是搜尋範圍,當 C t {C_t} Ct​越大時,接受B的機率越小,即目标點越不容易動。

是以, C t {C_t} Ct​是關于t的遞增函數,這樣就使前期盡可能地在全局搜尋,這樣就會防止陷入局部最優解;後期減小搜尋範圍,提高搜尋速度。(可以了解為渣男思想:前期廣撒網,後期精确突擊)

C t {C_t} Ct​的設定,模拟的是退火原理。

數學模組化暑期集訓23:模拟退火算法模拟退火算法類屬爬山法模拟退火算法流程 C t {C_t} Ct​的設定引入 C t {C_t} Ct​後的模拟退火算法流程例題1:求給定函數的最大值或者最小值例題2:旅行商問題總結
數學模組化暑期集訓23:模拟退火算法模拟退火算法類屬爬山法模拟退火算法流程 C t {C_t} Ct​的設定引入 C t {C_t} Ct​後的模拟退火算法流程例題1:求給定函數的最大值或者最小值例題2:旅行商問題總結
數學模組化暑期集訓23:模拟退火算法模拟退火算法類屬爬山法模拟退火算法流程 C t {C_t} Ct​的設定引入 C t {C_t} Ct​後的模拟退火算法流程例題1:求給定函數的最大值或者最小值例題2:旅行商問題總結

引入 C t {C_t} Ct​後的模拟退火算法流程

在定義好 C t {C_t} Ct​的表達式之後,可以将模拟退火算法的流程進行補充。

數學模組化暑期集訓23:模拟退火算法模拟退火算法類屬爬山法模拟退火算法流程 C t {C_t} Ct​的設定引入 C t {C_t} Ct​後的模拟退火算法流程例題1:求給定函數的最大值或者最小值例題2:旅行商問題總結

算法循環可以設定為2層,即第一次在高溫t的情況下進行周遊,之後逐漸降低溫度t,也就是改變 C t {C_t} Ct​,再次周遊。

例題1:求給定函數的最大值或者最小值

數學模組化暑期集訓23:模拟退火算法模拟退火算法類屬爬山法模拟退火算法流程 C t {C_t} Ct​的設定引入 C t {C_t} Ct​後的模拟退火算法流程例題1:求給定函數的最大值或者最小值例題2:旅行商問題總結

目标函數Obj_fun1

function y = Obj_fun1(x)
    y = 11*sin(x) + 7*cos(5*x);
end
           

模拟退火(含動畫示範)

%% SA 模拟退火: 求解函數y = 11*sin(x) + 7*cos(5*x)在[-3,3]内的最大值(動畫示範)
tic
clear; clc

%% 繪制函數的圖形
x = -3:0.1:3;
y = 11*sin(x) + 7*cos(5*x);
figure
plot(x,y,'b-')
hold on  % 不關閉圖形,繼續在上面畫圖

%% 參數初始化
narvs = 1; % 變量個數
T0 = 100;   % 初始溫度
T = T0; % 疊代中溫度會發生改變,第一次疊代時溫度就是T0
maxgen = 200;  % 最大疊代次數
Lk = 100;  % 每個溫度下的疊代次數
alfa = 0.95;  % 溫度衰減系數
x_lb = -3; % x的下界
x_ub = 3; % x的上界

%%  随機生成一個初始解
x0 = zeros(1,narvs);
for i = 1: narvs
    x0(i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(1);    
end
y0 = Obj_fun1(x0); % 計算目前解的函數值
h = scatter(x0,y0,'*r');  % scatter是繪制二維散點圖的函數(這裡傳回h是為了得到圖形的句柄,未來我們對其位置進行更新)

%% 定義一些儲存中間過程的量,友善輸出結果和畫圖
max_y = y0;     % 初始化找到的最佳的解對應的函數值為y0
MAXY = zeros(maxgen,1); % 記錄每一次外層循環結束後找到的max_y (友善畫圖)

%% 模拟退火過程
for iter = 1 : maxgen  % 外循環, 我這裡采用的是指定最大疊代次數
    for i = 1 : Lk  % 内循環,在每個溫度下開始疊代
        y = randn(1,narvs);  % 生成1行narvs列的N(0,1)随機數
        z = y / sqrt(sum(y.^2)); % 根據新解的産生規則計算z
        x_new = x0 + z*T; % 根據新解的産生規則計算x_new的值
        % 如果這個新解的位置超出了定義域,就對其進行調整
        for j = 1: narvs
            if x_new(j) < x_lb(j)
                r = rand(1);
                x_new(j) = r*x_lb(j)+(1-r)*x0(j);
            elseif x_new(j) > x_ub(j)
                r = rand(1);
                x_new(j) = r*x_ub(j)+(1-r)*x0(j);
            end
        end
        x1 = x_new;    % 将調整後的x_new指派給新解x1
        y1 = Obj_fun1(x1);  % 計算新解的函數值
        if y1 > y0    % 如果新解函數值大于目前解的函數值
            x0 = x1; % 更新目前解為新解
            y0 = y1;
        else
            p = exp(-(y0 - y1)/T); % 根據Metropolis準則計算一個機率
            if rand(1) < p   % 生成一個随機數和這個機率比較,如果該随機數小于這個機率
                x0 = x1; % 更新目前解為新解
                y0 = y1;
            end
        end
        % 判斷是否要更新找到的最佳的解
        if y0 > max_y  % 如果目前解更好,則對其進行更新
            max_y = y0;  % 更新最大的y
            best_x = x0;  % 更新找到的最好的x
        end
    end
    MAXY(iter) = max_y; % 儲存本輪外循環結束後找到的最大的y
    T = alfa*T;   % 溫度下降
    pause(0.01)  % 暫停一段時間(機關:秒)後再接着畫圖
    h.XData = x0;  % 更新散點圖句柄的x軸的資料(此時解的位置在圖上發生了變化)
    h.YData = Obj_fun1(x0); % 更新散點圖句柄的y軸的資料(此時解的位置在圖上發生了變化)
end

disp('最佳的位置是:'); disp(best_x)
disp('此時最優值是:'); disp(max_y)

pause(0.5)
h.XData = [];  h.YData = [];  % 将原來的散點删除
scatter(best_x,max_y,'*r');  % 在最大值處重新标上散點
title(['模拟退火找到的最大值為', num2str(max_y)])  % 加上圖的标題

%% 畫出每次疊代後找到的最大y的圖形
figure
plot(1:maxgen,MAXY,'b-');
xlabel('疊代次數');
ylabel('y的值');
toc
           

例題2:旅行商問題

數學模組化暑期集訓23:模拟退火算法模拟退火算法類屬爬山法模拟退火算法流程 C t {C_t} Ct​的設定引入 C t {C_t} Ct​後的模拟退火算法流程例題1:求給定函數的最大值或者最小值例題2:旅行商問題總結
tic
rng('shuffle')  % 控制随機數的生成,否則每次打開matlab得到的結果都一樣
% https://ww2.mathworks.cn/help/matlab/math/why-do-random-numbers-repeat-after-startup.html
% https://ww2.mathworks.cn/help/matlab/ref/rng.html
clear;clc
% 38個城市,TSP資料集網站(http://www.tsp.gatech.edu/world/djtour.html) 上公測的最優結果6656。
coord = [11003.611100,42102.500000;11108.611100,42373.888900;11133.333300,42885.833300;11155.833300,42712.500000;11183.333300,42933.333300;11297.500000,42853.333300;11310.277800,42929.444400;11416.666700,42983.333300;11423.888900,43000.277800;11438.333300,42057.222200;11461.111100,43252.777800;11485.555600,43187.222200;11503.055600,42855.277800;11511.388900,42106.388900;11522.222200,42841.944400;11569.444400,43136.666700;11583.333300,43150.000000;11595.000000,43148.055600;11600.000000,43150.000000;11690.555600,42686.666700;11715.833300,41836.111100;11751.111100,42814.444400;11770.277800,42651.944400;11785.277800,42884.444400;11822.777800,42673.611100;11846.944400,42660.555600;11963.055600,43290.555600;11973.055600,43026.111100;12058.333300,42195.555600;12149.444400,42477.500000;12286.944400,43355.555600;12300.000000,42433.333300;12355.833300,43156.388900;12363.333300,43189.166700;12372.777800,42711.388900;12386.666700,43334.722200;12421.666700,42895.555600;12645.000000,42973.333300];
n = size(coord,1);  % 城市的數目


figure  % 建立一個圖形視窗
plot(coord(:,1),coord(:,2),'o');   % 畫出城市的分布散點圖
hold on % 等一下要接着在這個圖形上畫圖的

d = zeros(n);   % 初始化兩個城市的距離矩陣全為0
for i = 2:n  
    for j = 1:i  
        coord_i = coord(i,:);   x_i = coord_i(1);     y_i = coord_i(2);  % 城市i的橫坐标為x_i,縱坐标為y_i
        coord_j = coord(j,:);   x_j = coord_j(1);     y_j = coord_j(2);  % 城市j的橫坐标為x_j,縱坐标為y_j
        d(i,j) = sqrt((x_i-x_j)^2 + (y_i-y_j)^2);   % 計算城市i和j的距離
    end
end
d = d+d';   % 生成距離矩陣的對稱的一面


%% 參數初始化
T0 = 1000;   % 初始溫度
T = T0; % 疊代中溫度會發生改變,第一次疊代時溫度就是T0
maxgen = 1000;  % 最大疊代次數
Lk = 500;  % 每個溫度下的疊代次數
alpfa = 0.95;  % 溫度衰減系數

%%  随機生成一個初始解
path0 = randperm(n);  % 生成一個1-n的随機打亂的序列作為初始的路徑
result0 = calculate_tsp_d(path0,d); % 調用我們自己寫的calculate_tsp_d函數計算目前路徑的距離

%% 定義一些儲存中間過程的量,友善輸出結果和畫圖
min_result = result0;     % 初始化找到的最佳的解對應的距離為result0
RESULT = zeros(maxgen,1); % 記錄每一次外層循環結束後找到的min_result (友善畫圖)

%% 模拟退火過程
for iter = 1 : maxgen  % 外循環, 我這裡采用的是指定最大疊代次數
    for i = 1 : Lk  %  内循環,在每個溫度下開始疊代
        path1 = gen_new_path(path0);  % 調用我們自己寫的gen_new_path函數生成新的路徑
        result1 = calculate_tsp_d(path1,d); % 計算新路徑的距離
        %如果新解距離短,則直接把新解的值指派給原來的解
        if result1 < result0    
            path0 = path1; % 更新目前路徑為新路徑
            result0 = result1; 
        else
            p = exp(-(result1 - result0)/T); % 根據Metropolis準則計算一個機率
            if rand(1) < p   % 生成一個随機數和這個機率比較,如果該随機數小于這個機率
                path0 = path1;  % 更新目前路徑為新路徑
                result0 = result1; 
            end
        end
        % 判斷是否要更新找到的最佳的解
        if result0 < min_result  % 如果目前解更好,則對其進行更新
            min_result = result0;  % 更新最小的距離
            best_path = path0;  % 更新找到的最短路徑
        end
    end
    RESULT(iter) = min_result; % 儲存本輪外循環結束後找到的最小距離
    T = alpfa*T;   % 溫度下降       
end


disp('最佳的方案是:'); disp(mat2str(best_path))
disp('此時最優值是:'); disp(min_result)


best_path = [best_path,best_path(1)];   % 在最短路徑的最後面加上一個元素,即第一個點(我們要生成一個封閉的圖形)
n = n+1;  % 城市的個數加一個(緊随着上一步)
for i = 1:n-1 
    j = i+1;
    coord_i = coord(best_path(i),:);   x_i = coord_i(1);     y_i = coord_i(2); 
    coord_j = coord(best_path(j),:);   x_j = coord_j(1);     y_j = coord_j(2);
    plot([x_i,x_j],[y_i,y_j],'-b')    % 每兩個點就作出一條線段,直到所有的城市都走完
%     pause(0.02)  % 暫停0.02s再畫下一條線段
    hold on
end

%% 畫出每次疊代後找到的最短路徑的圖形
figure
plot(1:maxgen,RESULT,'b-');
xlabel('疊代次數');
ylabel('最短路徑');

toc
           

總結

模拟退火作為智能優化算法(啟發式算法)中的一種,運作速度遠高于蒙特卡洛法。

缺點:matlab無通用内置函數,應對不同問題需要根據需要自行編寫代碼。

繼續閱讀