模擬退火算法 
　　模擬退火算法來源于固體退火原理，將固體加溫至充分高，再讓其徐徐冷卻，加溫時，固體內部粒子隨溫升變為無序狀，

內能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達到平衡態，最后在常溫時達到基態，內能減為最小。根據Metropolis

準則，粒子在溫度T時趨于平衡的概率為e-ΔE/(kT)，其中E為溫度T時的內能，ΔE為其改變量，k為Boltzmann常數。用固體退

火模擬組合優化問題，將內能E模擬為目標函數值f，溫度T演化成控制參數t，即得到解組合優化問題的模擬退火算法：由初始

解i和控制參數初值t開始，對當前解重復“產生新解→計算目標函數差→接受或舍棄”的迭代，并逐步衰減t值，算法終止時的

當前解即為所得近似最優解，這是基于蒙特卡羅迭代求解法的一種啟發式隨機搜索過程。退火過程由冷卻進度表(Cooling 

Schedule)控制，包括控制參數的初值t及其衰減因子Δt、每個t值時的迭代次數L和停止條件S。 
3.5.1 模擬退火算法的模型 
　　模擬退火算法可以分解為解空間、目標函數和初始解三部分。 
　模擬退火的基本思想: 
　　(1) 初始化：初始溫度T(充分大)，初始解狀態S(是算法迭代的起點)， 每個T值的迭代次數L 
　　(2) 對k=1，……，L做第(3)至第6步： 
　　(3) 產生新解S′ 
　　(4) 計算增量Δt′=C(S′)-C(S)，其中C(S)為評價函數 
　　(5) 若Δt′<0則接受S′作為新的當前解，否則以概率exp(-Δt′/T)接受S′作為新的當前解. 
　　(6) 如果滿足終止條件則輸出當前解作為最優解，結束程序。 
終止條件通常取為連續若干個新解都沒有被接受時終止算法。 
　　(7) T逐漸減少，且T->0，然后轉第2步。 
算法對應動態演示圖： 
模擬退火算法新解的產生和接受可分為如下四個步驟： 
　　第一步是由一個產生函數從當前解產生一個位于解空間的新解；為便于后續的計算和接受，減少算法耗時，通常選擇由當

前新解經過簡單地變換即可產生新解的方法，如對構成新解的全部或部分元素進行置換、互換等，注意到產生新解的變換方法

決定了當前新解的鄰域結構，因而對冷卻進度表的選取有一定的影響。 
　　第二步是計算與新解所對應的目標函數差。因為目標函數差僅由變換部分產生，所以目標函數差的計算最好按增量計算。

事實表明，對大多數應用而言，這是計算目標函數差的最快方法。 
　　第三步是判斷新解是否被接受,判斷的依據是一個接受準則，最常用的接受準則是Metropo1is準則: 若Δt′<0則接受S′作

為新的當前解S，否則以概率exp(-Δt′/T)接受S′作為新的當前解S。 
　　第四步是當新解被確定接受時，用新解代替當前解，這只需將當前解中對應于產生新解時的變換部分予以實現，同時修正

目標函數值即可。此時，當前解實現了一次迭代。可在此基礎上開始下一輪試驗。而當新解被判定為舍棄時，則在原當前解的

基礎上繼續下一輪試驗。 
　　模擬退火算法與初始值無關，算法求得的解與初始解狀態S(是算法迭代的起點)無關；模擬退火算法具有漸近收斂性，已在

理論上被證明是一種以概率l 收斂于全局最優解的全局優化算法；模擬退火算法具有并行性。 

3.5.2 模擬退火算法的簡單應用 
　　作為模擬退火算法應用，討論貨郎擔問題(Travelling Salesman Problem，簡記為TSP)：設有n個城市，用數碼1,…,n代表

。城市i和城市j之間的距離為d(i，j) i, j=1,…,n．TSP問題是要找遍訪每個域市恰好一次的一條回路，且其路徑總長度為最

短.。 
　　求解TSP的模擬退火算法模型可描述如下： 
　　解空間 解空間S是遍訪每個城市恰好一次的所有回路，是{1，……，n}的所有循環排列的集合，S中的成員記為(w1,w2 ,…

…，wn)，并記wn+1= w1。初始解可選為(1，……，n) 
　　目標函數 此時的目標函數即為訪問所有城市的路徑總長度或稱為代價函數： 

　　我們要求此代價函數的最小值。 
　　新解的產生 隨機產生1和n之間的兩相異數k和m，若k<m，則將 
　　(w1, w2 ,…，wk , wk+1 ,…，wm ,…，wn) 
　　變為： 
　　(w1, w2 ,…，wm , wm-1 ,…，wk+1 , wk ,…，wn). 
　　如果是k>m，則將 
　　(w1, w2 ,…，wk , wk+1 ,…，wm ,…，wn) 
　　變為： 
　　(wm, wm-1 ,…，w1 , wm+1 ,…，wk-1 ,wn , wn-1 ,…，wk). 
　　上述變換方法可簡單說成是“逆轉中間或者逆轉兩端”。 
　　也可以采用其他的變換方法，有些變換有獨特的優越性，有時也將它們交替使用，得到一種更好方法。 
　　代價函數差 設將(w1, w2 ,……，wn)變換為(u1, u2 ,……，un), 則代價函數差為： 

根據上述分析，可寫出用模擬退火算法求解TSP問題的偽程序： 
Procedure TSPSA: 
　begin 
　　init-of-T; { T為初始溫度} 
　　S={1，……，n}; {S為初始值} 
　　termination=false; 
　　while termination=false 
　　　begin 
　　　　for i=1 to L do 
　　　　　　begin 
　　　　　　　　generate(S′form S); { 從當前回路S產生新回路S′} 
　　　　　　　　Δt:=f(S′))-f(S);{f(S)為路徑總長} 
　　　　　　　　IF(Δt<0) OR (EXP(-Δt/T)>Random-of-[0,1]) 
　　　　　　　　S=S′; 
　　　　　　　　IF the-halt-condition-is-TRUE THEN 
　　　　　　　　termination=true; 
　　　　　　End; 
　　　　T_lower; 
　　　End; 
　End 
　　模擬退火算法的應用很廣泛，可以較高的效率求解最大截問題(Max Cut Problem)、0-1背包問題(Zero One Knapsack 

Problem)、圖著色問題(Graph Colouring Problem)、調度問題(Scheduling Problem)等等。 

3.5.3 模擬退火算法的參數控制問題 
　　模擬退火算法的應用很廣泛，可以求解NP完全問題，但其參數難以控制，其主要問題有以下三點： 
　　(1) 溫度T的初始值設置問題。 
　　溫度T的初始值設置是影響模擬退火算法全局搜索性能的重要因素之一、初始溫度高，則搜索到全局最優解的可能性大，但

因此要花費大量的計算時間；反之，則可節約計算時間，但全局搜索性能可能受到影響。實際應用過程中，初始溫度一般需要

依據實驗結果進行若干次調整。 
　　(2) 退火速度問題。 
　　模擬退火算法的全局搜索性能也與退火速度密切相關。一般來說，同一溫度下的“充分”搜索(退火)是相當必要的，但這

需要計算時間。實際應用中，要針對具體問題的性質和特征設置合理的退火平衡條件。 
　　(3) 溫度管理問題。 
　　溫度管理問題也是模擬退火算法難以處理的問題之一。實際應用中，由于必須考慮計算復雜度的切實可行性等問題，常采

用如下所示的降溫方式： 

T(t+1)＝k×T(t) 
式中k為正的略小于1.00的常數，t為降溫的次數 

使用SA解決TSP問題的Matlab程序：

function out = tsp(loc)
% TSP Traveling salesman problem (TSP) using SA (simulated annealing).
% TSP by itself will generate 20 cities within a unit cube and
% then use SA to slove this problem.
%
% TSP(LOC) solve the traveling salesman problem with cities'
% coordinates given by LOC, which is an M by 2 matrix and M is
% the number of cities.
%
% For example:
%
% loc = rand(50, 2);
% tsp(loc);
if nargin == 0,
% The following data is from the post by Jennifer Myers (jmyers@nwu.
edu)
edu)
% to comp.ai.neural-nets. It's obtained from the figure in
% Hopfield & Tank's 1985 paper in Biological Cybernetics
% (Vol 52, pp. 141-152).
loc = [0.3663, 0.9076; 0.7459, 0.8713; 0.4521, 0.8465;
0.7624, 0.7459; 0.7096, 0.7228; 0.0710, 0.7426;
0.4224, 0.7129; 0.5908, 0.6931; 0.3201, 0.6403;
0.5974, 0.6436; 0.3630, 0.5908; 0.6700, 0.5908;
0.6172, 0.5495; 0.6667, 0.5446; 0.1980, 0.4686;
0.3498, 0.4488; 0.2673, 0.4274; 0.9439, 0.4208;
0.8218, 0.3795; 0.3729, 0.2690; 0.6073, 0.2640;
0.4158, 0.2475; 0.5990, 0.2261; 0.3927, 0.1947;
0.5347, 0.1898; 0.3960, 0.1320; 0.6287, 0.0842;
0.5000, 0.0396; 0.9802, 0.0182; 0.6832, 0.8515];
end
NumCity = length(loc); % Number of cities
distance = zeros(NumCity); % Initialize a distance matrix
% Fill the distance matrix
for i = 1:NumCity,
for j = 1:NumCity,
distance(i, j) = norm(loc(i, - loc(j, );
distance(i, j) = norm(loc(i, - loc(j, );
end
end
% To generate energy (objective function) from path
%path = randperm(NumCity);
%energy = sum(distance((path-1)*NumCity + [path(2:NumCity) path(1)]));
% Find typical values of dE
count = 20;
all_dE = zeros(count, 1);
for i = 1:count
path = randperm(NumCity);
energy = sum(distance((path-1)*NumCity + [path(2:NumCity)
path(1)]));
new_path = path;
index = round(rand(2,1)*NumCity+.5);
inversion_index = (min(index):max(index));
new_path(inversion_index) = fliplr(path(inversion_index));
all_dE(i) = abs(energy - ...
sum(sum(diff(loc([new_path new_path(1)],)'.^2)));
end
dE = max(all_dE);
dE = max(all_dE);
temp = 10*dE; % Choose the temperature to be large enough
fprintf('Initial energy = %f\n\n',energy);
% Initial plots
out = [path path(1)];
plot(loc(out(, 1), loc(out(, 2),'r.', 'Markersize', 20);
axis square; hold on
h = plot(loc(out(, 1), loc(out(, 2)); hold off
MaxTrialN = NumCity*100; % Max. # of trials at a
temperature
MaxAcceptN = NumCity*10; % Max. # of acceptances at a
temperature
StopTolerance = 0.005; % Stopping tolerance
TempRatio = 0.5; % Temperature decrease ratio
minE = inf; % Initial value for min. energy
maxE = -1; % Initial value for max. energy
% Major annealing loop
while (maxE - minE)/maxE > StopTolerance,
minE = inf;
minE = inf;
maxE = 0;
TrialN = 0; % Number of trial moves
AcceptN = 0; % Number of actual moves
while TrialN < MaxTrialN & AcceptN < MaxAcceptN,
new_path = path;
index = round(rand(2,1)*NumCity+.5);
inversion_index = (min(index):max(index));
new_path(inversion_index) =
fliplr(path(inversion_index));
new_energy = sum(distance( ...
(new_path-1)*NumCity+[new_path(2:NumCity)
new_path(1)]));
if rand < exp((energy - new_energy)/temp), % 
accept
it!
energy = new_energy;
path = new_path;
minE = min(minE, energy);
maxE = max(maxE, energy);
AcceptN = AcceptN + 1;
end
TrialN = TrialN + 1;
end
end
% Update plot
out = [path path(1)];
set(h, 'xdata', loc(out(, 1), 'ydata', loc(out(, 2));
drawnow;
% Print information in command window
fprintf('temp. = %f\n', temp);
tmp = sprintf('%d ',path);
fprintf('path = %s\n', tmp);
fprintf('energy = %f\n', energy);
fprintf('[minE maxE] = [%f %f]\n', minE, maxE);
fprintf('[AcceptN TrialN] = [%d %d]\n\n', AcceptN, TrialN);
% Lower the temperature
temp = temp*TempRatio;
end
% Print sequential numbers in the graphic window
for i = 1:NumCity,
text(loc(path(i),1)+0.01, loc(path(i),2)+0.01, int2str(i), ...
'fontsize', 8);
end