模擬退火演算法學習筆記

最開始聽說到退火這個概念，是在高中的時候很喜歡看的那些量子物理相關的文章里。

量子物理最初給我留下的印象是非常神秘深奧和捉摸不透的，但是其實只是因為不了解。而模擬退火演算法也是這樣。剛聽到這個名字，會覺得，哇丟，賊高大上。不過學習之後，會發現這是一個非常有趣，而且淺顯易懂的東西。

這篇文章，我主要用 4 個不同類型的例子，使用 MATLAB，來寫下我自己學習模擬退火演算法之後的理解。

模擬退火演算法主要用於難以準確求出具體的解的問題之中。通過多次迭代，它可以不斷地接近最優解。

只是接近最優解？

是的！對於一些問題，如TSP問題，如果把所有可能的解都遍歷一遍，需要的時間是隨著城市的數量增加而呈爆炸性增長的，因此對於多個城市的TSP問題，老老實實去解，會花費特別多的時間。若是對於解的精度沒有太高的要求，則可以損失一些精度來換取速度。

舉個簡單的例子，比如我們要找到方程 的最小值，雖然對於我們來說一眼就可以看出來解是 ，但是對於計算機呢？如果是更複雜的方程呢？

模擬退火演算法的核心思想是：首先隨機選擇一個解作為開始，接下來產生一個隨機擾動，如果找到比上一個解更接近最優解的解，那麼就直接接受這個解。而如果找到的解離得更遠了，沒關係，以一定的概率接受。

這段文字可以說是很繞了，那奏拿上面的方程作為例子，一步一步拆解開，先看看會是什麼樣子叭！

首先我們限定一個區間，比如 ，然後在這個區間里，隨機選擇一個點， ，以及它對應的函數的數值 ，然後以這個點作為起點。

接下來產生一個隨機擾動，擾動的大小可以自己設計一些參數來調節。 ，並計算出該點對應的函數值 。

對兩個數值求差值，就可以判斷新的解是否更接近最優解。

如果dE小於0，則說明新的接更加接近我們想要的解，即接受該解。

我們會想，那更接近的話，接受，不就行了唄？為什麼還需要在離最優解更遠的時候也一定的概率接受呢？

如果函數不是 ，而是 [1]呢，哈哈哈哈開個玩笑。意思就是如果遇到一個有多個極值（或者說局部最優解）的函數，那麼如果不接受一個離最優解更遠的解，如果不是因為一開始隨機的點就落在了最優解所在的區間里，那麼就永遠也沒辦法得到最優解了。 以一定的概率去接受一個比較差勁的解的目的就是，讓它有機會跳出局部最優解而奔向全局最優解。

這個概率一般取 ，其中dE為當前解與上一個解的差值，k為常數，T為當前溫度。

e的負指數次冪小於1，因此可以十分方便地直接將p和範圍為(0, 1)的隨機數做比較，若p大於隨機數則接受該解。隨著溫度的降低，p的值也逐漸降低，於是接受更差的解的概率也就越小，退火過程就趨於穩定。

模擬退火的基本思想就這麼多了，接下來通過實例來理解這個過程8！

首先是模擬退火的精髓，判斷函數。

function [y] = judge(dE, t) if(dE < 0) y = 1; else d = exp(-(dE / t)); if(d > rand) y = 1; else y = 0; end endend

函數傳入兩個參數，dE與t，然後返回是否接受該解。

1、使用模擬退火演算法求解[0， 9]範圍內函數的最大值。

這裡就拿上面說到的那個又是cos又是sin的函數吧哈哈哈哈哈

clear;%範圍section_l = 0;section_h = 9;%繪製函數圖像draw_base();%初始溫度，停止溫度與降溫係數tmp = 1e5;tmp_min = 1e-3;alpha = 0.98;%生成初始隨機解x_old = (section_h - section_l) * rand() + section_l;x_new = x_old;s_old = val(x_old);s_new = s_old;text_lt = text(0, 0, Init);%計數器counter = 0;%退火的主體部分，一個循環while(tmp > tmp_min) %隨機擾動 delta = (rand() - 0.5) * 3; x_new = x_old + delta; %擾動的值小於一半的區間範圍時，可以用這種辦法防止新值超出區間範圍 if(x_new < section_l || x_new > section_h) x_new = x_new - 2 * delta; end s_new = val(x_new); %求差值，這裡是找最大值而非最小值，所以不是s_new - s_old dE = s_old - s_new; %判斷 j = judge(dE, tmp); if(j) s_old = s_new; x_old = x_new; end %只有當dE < 0的情況下才降溫 if(dE < 0) delete(text_lt); hold on, scatter(x_old, s_old); hold on, text_lt = text(0.3, 21, [tmp: , num2str(tmp)]); pause(0.01); %上面是繪圖的代碼，下面降溫 tmp = tmp * alpha; else counter = counter + 1; end %當接受更差的解的概率太小時，若又長時間找不到更優的解，那麼退出循環，結束演算法 if(counter > 10000) break; end endfunction [y] = val(x) y = x + 10 * sin(5 * x) + 7 * cos(4 * x);endfunction draw_base() X = 0: 0.01:9; M = val(X); plot(X, M);end

整個程序可以大致分為幾個部分：

函數的值，val()函數，傳入自變數，輸出函數值。判斷函數，輸入dE與T，輸出是否接受新解。最後是退火的主體部分，在一個while循環中，首先設定退出循環的條件，比如在本例中是退火達到了某個最低的溫度，或者一定次數下都只能找到更劣質的解。

2、使用模擬退火演算法求解到點集直線距離之和最短的點

clear;global num_of_dots;global dots_x;global dots_y;num_of_dots = 12;dots_x = [0 1.5 1.5 -2 -4 -5 2 4 5 -2 -3 3 -1.5 -2.5 0.1];dots_y = [0 10 12 11 -8 2 -1.5 -2.5 1 -2 8 6 0 0 0.2];draw_dist(), hold on;scatter(dots_x, dots_y, MarkerFaceColor,blue);x_min = min(dots_x);y_min = min(dots_y);x_max = max(dots_x);y_max = max(dots_y);x_range = x_max - x_min;y_range = y_max - y_min;tmp = 1000;tmp_max = tmp;tmp_min = 0.001;alpha = 0.98;x_old = rand() * x_range + x_min;y_old = rand() * y_range + y_min;x_new = x_old;y_new = y_old;s_old = dist(x_old, y_old);s_new = s_old;txttt = text(0, 0, Init);while(tmp > tmp_min) %因為是在二維平面里，所以產生兩個維度的隨機擾動 delta1 = (rand - 0.5); delta2 = (rand - 0.5); x_new = delta1 + x_old; y_new = delta2 + y_old; if(x_new > x_max || x_new < x_min) x_new = x_new - 2 * delta1; end if(y_new > y_max || y_new < y_min) y_new = y_new - 2 * delta2; end s_new = dist(x_new, y_new); dE = s_new - s_old; j = judge(dE, tmp); if(j) s_old = s_new; x_old = x_new; y_old = y_new; end delete(txttt); txttt = text(-4, 6, {[Dist: ,num2str(s_old)];[Tmp: , num2str(tmp)]}); hold on, scatter(x_old, y_old, .), pause(0.01); tmp = tmp * alpha;endtext(x_old, y_old, [Solve: , num2str(s_old)]);%計算點到點集的距離function [d] = dist(x, y) global num_of_dots; global dots_x; global dots_y; d = 0; for ii = 1:num_of_dots d = d + sqrt((x - dots_x(ii))^2 + (y - dots_y(ii))^2); endendfunction draw_dist() global dots_x; global dots_y; x_min = floor(min(dots_x)); y_min = floor(min(dots_y)); x_max = ceil(max(dots_x)); y_max = ceil(max(dots_y)); resol = 10; dx = x_min: 1/resol : x_max; dy = y_min: 1/resol : y_max; [xx, yy] = meshgrid(dx, dy); zz = zeros(length(dy), length(dx)); for ii = 1:length(dx) for jj = 1:length(dy) zz(jj, ii) = dist(dx(ii), dy(jj)); end end contour(xx, yy, zz);end

這個問題其實使用爬山演算法，梯度下降等也可以很方便地求解，因為它只有一個極值

3、TSP問題[2]

賊炫酷有沒有！

clear;cities1=[0.9695,0.6606,0.5906,0.2124,0.0398,0.1367,0.9536,0.6091,0.8767,0.8148,0.3876,0.7041,0.0213,0.3429,0.7471,0.4606,0.7695,0.5006,0.3124,0.0098,0.3637,0.5336,0.2091,0.4767,0.4148,0.5876,0.6041,0.3213,0.6429,0.7471; 0.6740,0.9500,0.5029,0.8274,0.9697,0.5979,0.2184,0.7148,0.2395,0.2867,0.8200,0.3296,0.1649,0.3025,0.8192,0.6500,0.7420,0.0229,0.7274,0.4697,0.0979,0.2684,0.7948,0.4395,0.8867,0.3200,0.5296,0.3649,0.7025,0.9192];tmp = 1e3;tmp_min = 1e-5;alpha = 0.99;s_old = dist(cities1);cit_new = cities1;s_new = s_old;txttt = text(0, 0, Init);counter = 0;while(tmp > tmp_min) %產生隨機擾動，對於TSP問題，隨機擾動一般為交換幾座城市的順序 cit_new = new_solve(cities1); s_new = dist(cit_new); dE = (s_new - s_old) * 500; j = judge(dE, tmp); if(j) cities1 = cit_new; s_old = s_new; end if(dE < 0) counter = 0; %降溫 tmp = tmp * alpha; %繪圖 clf; text(0.8, 0.9, {[tmp: , num2str(tmp)];[dist: , num2str(s_old)]}); hold on; draw_route(cities1); pause(0.01); else counter =counter + 1; end if(counter > 1e4) break; endendfunction draw_route(in) plot(in(1,:), in(2,:),-o);endfunction [d] = dist(in) num = size(in, 2) - 1; d = 0; for ii = 1: num d = d + sqrt(sum((in(:, ii+1) - in(:, ii)).^2)); endendfunction [cities2] = new_solve(in) cities2 = in; num_of_cities = size(cities2,2); while 1 cit1 = ceil(rand() * num_of_cities); cit2 = ceil(rand() * num_of_cities); if(cit1 ~= cit2) break; end end tmp = cities2(:, cit1); cities2(:, cit1) = cities2(:, cit2); cities2(:, cit2) = tmp;end

4、應用題

clear;v_old = rand();tao_old = rand();v_new = v_old;tao_new = tao_old;s_old = vt(v_old, tao_old);s_new = s_old;tmp = 1e5;tmp_min = 1e-3;alpha = 0.98;draw_vt();txttt = text(0, 0, Init);while(tmp > tmp_min) delta1 = (rand - 0.5); delta2 = (rand - 0.5); v_new = delta1 + v_old; tao_new = delta2 + tao_old; if(tao_new == 0) tao_new = tao_new + delta2; end s_new = vt(v_new, tao_new); dE = s_new - s_old; j = judge(dE, tmp); if(j) s_old = s_new; v_old = v_new; tao_old = tao_new; end delete(txttt); hold on, scatter(v_old, tao_old, .); txttt = text(-20, 2.5, {[diff: ,num2str(s_old)];[Tmp: , num2str(tmp)];[Tao: , num2str(tao_old)];[V0: , num2str(v_old)]},Color, white); pause(0.01); tmp = tmp * alpha;endfunction draw_vt() V0 = linspace(-30, 30, 50); Tao = linspace(0.1, 5, 50); zz = zeros(length(V0)); for ii = 1:length(V0) for jj = 1:length(Tao) zz(jj, ii) = abs(6.1434 - vt(V0(ii), Tao(jj))); end end gca = pcolor(V0, Tao, zz); set(gca, LineStyle,none); shading interp;endfunction [y] = vt(V0, Tao) V = 14; T = [0.3 0.5 1.0 2.0 4.0 7.0]; res_origin = [5.6873 6.1434 7.1633 8.8626 11.0328 12.6962]; res = V - (V - V0) * exp(-T/Tao); delta = res - res_origin; y = sum(delta.*delta);end

[1] 這個函數來自 如何通俗易懂地解釋遺傳演算法？有什麼例子？

[2] Aravind Seshadri.Simulated Annealing For Traveling Salesman Problem