前言
模擬退火算法（SA）是較為常見的現代優化算法之一，常用於旅行商（TSP）問題中。數學建模里學生們常常使用該算法，甚至是為了使用這個算法而使用這個算法，讓評委老師們審美疲勞。評委老師明確表明使用所謂"神算法"（神經網絡，模擬退火，遺傳算法等等）而過於牽強者拿不了高分（見：http://special.univs.cn/service/jianmo/sxjmyw/2018/1128/1187951_15.shtml）。希望大家不要覺得它名詞高級就認為它能吸引評委眼睛，評委畢竟是教授，不可能被幾個名詞唬住。
但是呢，我們是學生，不能因為它不能隨便用就不學習它，而在編程的環節中，我們亦有收獲，況且愛因斯坦也是從一加一開始學起的，所以模擬退火算法還是有學習的必要的。話說的有點多，下面進入主題。
算法框架
模擬退火算法可以粗分為以下幾個步驟：
1，初始溫度的設置、初始解的生成、設置每個溫度下產生解的個數。
2，產生新解。
3，計算代價函數差。
4，Metropolis判別。（別被名詞嚇住，形式上是很簡單的一個原則）
5，降溫。
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.