為了參加國賽，這幾天學了模擬退火演算法，整理下當做模板方便國賽的時候用。

模擬退火用於處理最優化問題，可以求出當目標函式取得最小值時的決策變數的值。

在編寫程式時需要根據具體問題設計演算法，演算法描述為：

（1）解空間（初始解）

（2）目標函式

（3）新解的產生
          ① 2 變換法

          ② 3 變換法

（4）代價函式差

（5）接受準則

（6）降溫

（7）結束條件

下面MATLAB程式用於求解非線性規劃：

min f(x)=x1^2+x2^2+8

st.

x1^2-x2>=0

-x1-x2^2+2=0

x1,x2>=0
 

MATLAB程式碼：

clear
clc

%生成初始解
sol_new2=1;%（1）解空間（初始解）
sol_new1=2-sol_new2^2;
sol_current1 = sol_new1; 
sol_best1 = sol_new1;
sol_current2 = sol_new2; 
sol_best2 = sol_new2;
E_current = inf;
E_best = inf;

rand('state',sum(clock)); %初始化隨機數發生器
t=90; %初始溫度
tf=89.9; %結束溫度
a = 0.99; %溫度下降比例

while t>=tf%（7）結束條件
    for r=1:1000 %退火次數
        
        %產生隨機擾動（3）新解的產生
        sol_new2=sol_new2+rand*0.2;
        sol_new1=2-sol_new2^2;
        
        %檢查是否滿足約束
        if sol_new1^2-sol_new2>=0 && -sol_new1-sol_new2^2+2==0 && sol_new1>=0 &&sol_new2>=0
        else
            sol_new2=rand*2;
            sol_new1=2-sol_new2^2;
            continue;
        end
        
        %退火過程
        E_new=sol_new1^2+sol_new2^2+8;%（2）目標函式
        if E_new<E_current%（5）接受準則
                E_current=E_new;
                sol_current1=sol_new1;
                sol_current2=sol_new2;
                if E_new<E_best
                    %把冷卻過程中最好的解儲存下來
                    E_best=E_new;
                    sol_best1=sol_new1;
                    sol_best2=sol_new2;
                end
        else
                if rand<exp(-(E_new-E_current)/t)%（4）代價函式差
                    E_current=E_new;
                    sol_current1=sol_new1;
                    sol_current2=sol_new2;
                else
                    sol_new1=sol_current1;
                    sol_new2=sol_current2;
                end
        end
        plot(r,E_best,'*')
        hold on
    end
    t=t*a;%（6）降溫
end

disp('最優解為：')
disp(sol_best1)
disp(sol_best2)
disp('目標表達式的最小值等於：')
disp(E_best)
 

執行結果：

最優解為：
    1.0038

    0.9981

目標表達式的最小值等於：
   10.0038

司守奎的演算法大全給出了利用模擬退火求解TSP問題的MATLAB演算法：

clc,clear
load sj.txt %載入敵方100 個目標的資料，資料按照表格中的位置儲存在純文字檔案sj.txt 中
x=sj(:,1:2:8);x=x(:);
y=sj(:,2:2:8);y=y(:);
sj=[x y];
d1=[70,40];
sj=[d1;sj;d1];
sj=sj*pi/180;
%距離矩陣d
d=zeros(102);
for i=1:101
    for j=i+1:102
        temp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2));
        d(i,j)=6370*acos(temp);
    end
end
d=d+d';
S0=[];Sum=inf;
rand('state',sum(clock));
for j=1:1000
    S=[1 1+randperm(100),102];
    temp=0;
    for i=1:101
        temp=temp+d(S(i),S(i+1));
    end
    if temp<Sum
        S0=S;Sum=temp;
    end
end
e=0.1^30;L=20000;at=0.999;T=1;
%退火過程
for k=1:L
    %產生新解
    c=2+floor(100*rand(1,2));
    c=sort(c);
    c1=c(1);c2=c(2);
    %計算代價函式值
    df=d(S0(c1-1),S0(c2))+d(S0(c1),S0(c2+1))-d(S0(c1-1),S0(c1))-        d(S0(c2),S0(c2+1));
    %接受準則
    if df<0
        S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
        Sum=Sum+df;
    elseif exp(-df/T)>rand(1)
        S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
        Sum=Sum+df;
    end
    T=T*at;
    if T<e
        break;
    end
end
% 輸出巡航路徑及路徑長度
S0,Sum