前言html
模擬退火算法（SA）是較爲常見的現代優化算法之一，經常使用於旅行商（TSP）問題中。數學建模裏學生們經常使用該算法，甚至是爲了使用這個算法而使用這個算法，讓評委老師們審美疲勞。評委老師明確代表使用所謂"神算法"（神經網絡，模擬退火，遺傳算法等等）而過於牽強者拿不了高分（見：http://special.univs.cn/service/jianmo/sxjmyw/2018/1128/1187951_15.shtml）。但願你們不要以爲它名詞高級就認爲它能吸引評委眼睛，評委畢竟是教授，不可能被幾個名詞唬住。算法
可是呢，咱們是學生，不能由於它不能隨便用就不學習它，而在編程的環節中，咱們亦有收穫，何況愛因斯坦也是從一加一開始學起的，因此模擬退火算法仍是有學習的必要的。話說的有點多，下面進入主題。編程
算法框架網絡
模擬退火算法能夠粗分爲如下幾個步驟：框架
1，初始溫度的設置、初始解的生成、設置每一個溫度下產生解的個數。函數
2，產生新解。學習
3，計算代價函數差。測試
4，Metropolis判別。（別被名詞嚇住，形式上是很簡單的一個原則）優化
5，降溫。spa
6，判斷溫度是否小於一個給定量。是，則結束；否，則跳轉到第2步。
如下對每一個步驟作詳細的解釋。
初始溫度的設置、初始解的生成、設置每一個溫度下產生解的個數：
初始溫度與降溫係數、終止溫度息息相關，它們仨決定了迭代的次數，具體公式爲：終止溫度 ≈ 初始溫度*(降溫係數)^迭代次數，其中1降溫係數是(0,1)的一個常數 。對於小規模的問題，通常初始溫度取1000，終止溫度取1e-4，降溫係數取0.9。這些參數直接影響了Metropolis判別。
初始解的生成通常是隨機生成，由於SA算法對初始解的依賴並不大。
設置每一個溫度下產生解的個數的目的是增長樣本量，加強算法的健壯性。多組解稱爲Metropolis的鏈長。
產生新解
產生新解的方式有不少種，最簡單的是隨機產生新解，固然也能夠按照必定的算法產生新解。隨機產生新解的好處是容易尋找到全局最優勢，按必定算法產生新解的好處是在知道解大體位置的前提下能夠更好的收斂。
計算代價函數差
計算代價函數差，即將最新解、最優解作差，獲得增量df。
Metropolis判別
該判別的公式以下：
意義爲，當代價函數差小於0（即新解比最優解的值還要小），接受最新解變爲最優解；當代價函數不小於0時，以必定機率P接受該解做爲最優解。由此能夠看到，溫度越高，劣解接受率也越低，換言之就是全局搜索能力更強。因此設置初始溫度時，要將溫度設置的比較高，避免陷入局部最優。
降溫
降溫方式較爲簡單，即新溫度=當前溫度*降溫係數。從中能夠看到迭代次數與初始溫度、降溫係數都有關係，但注意，初始溫度影響了Metropolis判別，因此若要想提升迭代次數，提升降溫係數較爲合適。最後就是判斷當前溫度是否小於終止溫度。以上就是關於模擬退火算法流程的我的理解，接下來將設計程序，以解決TSP問題。
程序設計
%模擬退火算法
%TSP問題
clc,clear
T0 =500;%初始溫度
Tend = 1e-4; %終止溫度
L=200; %各個溫度下的鏈長
q=0.98; %降溫速率
cityNum = 30; %城市個數
City = rand(cityNum,2); %隨機生成城市的座標
%%%%%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%%
D = Distance(City);
N = cityNum;
S1 = randperm(N); %產生一個隨機解
%解方程，計算迭代次數
Time = ceil(double(solve(['1000*(0.9)^x',num2str(Tend)])));
count = 0; %迭代計數器
Obj = zeros(Time,1); %每代路徑和
track = zeros(Time,N); %每代最優解
%迭代
while T0>Tend
count = count + 1;
temp = zeros(L,N+1);
for k = 1:L
%產生L組新解
S2 = newSolution(S1);
[S1,R] = Metropolis(S1,S2,D,T0);
temp(k,:)=[S1 R];
end
%記錄每次迭代過程的最優路線
[d0,index] = min(temp(:,end));
if count == 1 || d0 <Obj(count-1)
Obj(count) = d0;
else
Obj(count)=Obj(count - 1);
end
track(count,:)=temp(index,1:end-1);
T0 = q*T0;
end
fprintf('迭代次數：%d\n',count);
fprintf('最短路徑：%f\n',Obj(end));
%迭代圖
figure
plot(1:count,Obj);
xlabel('迭代次數');
ylabel('距離');
title('優化過程');
DrawPath(track(end,:),City);
%%%%%%%%%%%%%%%%%子函數，共5個%%%%%%%%%%%%%%%%%%%%%
%函數功能：求解距離矩陣
%輸入：城市座標 輸出：距離矩陣
function D = Distance(a)
r = size(a,1);
D = zeros(r,r);
for i = 1:r
for j = i+1:r
D(i,j)=sqrt((a(i,1)-a(j,1))^2 + (a(i,2)-a(j,2))^2);
D(j,i) = D(i,j);
end
end
end
%函數功能：生成新解
%輸入：舊解 輸出：新解
%隨機選兩個位置交換。
function S2 = newSolution(S1)
N = length(S1);
S2 = S1;
a = round( rand(1,2)*(N-1)+1); %產生兩個隨機位置用來交換
temp = S2(a(1));
S2(a(1))=S2(a(2));
S2(a(2)) = temp;
end
%Metropolis準則函數
%輸入：新解，舊解，距離矩陣，當前溫度
%輸出：判斷後的解，解的路徑
function [S,R] = Metropolis(S1,S2,D,T)
R1 = pathLength(D,S1);
R2 = pathLength(D,S2);
dT = R2 - R1;
if dT < 0
S=S2;
R=R2;
elseif exp(-dT/T) >= rand
S=S2;
R=R2;
else
S=S1;
R=R1;
end
end
%畫出路線軌跡圖
%輸入：最優路徑，城市座標 輸出：路徑圖
function DrawPath(BestTrack,City)
N = size(BestTrack,2);
cityArray = zeros(N,2);
for i = 1:N
cityArray(i,1) = City(BestTrack(i),1);
cityArray(i,2) = City(BestTrack(i),2);
end
figure;
hold on
plot(cityArray(:,1),cityArray(:,2),'-o','color',[0.5 0.5 0.5]);
plot([cityArray(1,1),cityArray(end,1)],[cityArray(1,2),cityArray(end,2)],...
'-','color',[0.5 0.5 0.5]);
hold off
xlabel('橫座標');
ylabel('縱座標');
title('TSP圖');
box on
end
%計算路線長度的函數
%輸入：距離矩陣，行走的順序 輸出：路徑和
function len = pathLength(D,S)
[r,c] = size(D);
NIND = size(S,1);
for i = 1:NIND
p = [S(i,:) S(i,1)];
i1 = p(1:end-1);
i2 = p(2:end);
len(i,1) = sum(D((i1-1)*c + i2));
end
end
注意子函數newSolution，個人策略是隨機挑選兩個位置進行交換。事實上還能夠採用其餘策略來產生新解，但我太懶了，因此沒想更好的產生新解的方法。
程序結果
最後，我測試了一下當城市數爲50，須要增大迭代次數，增大初始溫度才能獲得較好的解。固然，我產生新解的方式也較爲簡單，也能夠從這方面考慮，用更好的想法來產生新解。
參考文獻
[1] 鬱磊，史峯，王輝等，《matlab智能算法30個案例分析（第二版）》，2015.8.