模擬退火算法

模擬退火（SA）

物理過程由以下三個部分組成

1.加溫過程 問題的初始解

2.等溫過程 對應算法的Metropolis抽樣的過程

3.冷卻過程 控制參數的下降

默認的模擬退火是一個求最小值的過程，其中Metropolis准則是SA算法收斂於全局最優解的關鍵所在，Metropolis准則以一定的概率接受惡化解，這樣就使算法跳離局部最優的陷進

 

 1.模擬退火算法求解一元函數最值問題

使用simulannealbnd - Simulated annealing algorithm工具箱

求y=sin(10*pi*x)./x;在[1,2]的最值

下圖是用畫圖法求出最值的

x=1:0.01:2;
y=sin(10*pi*x)./x;
figure
hold on
plot(x,y,'linewidth',1.5);
ylim([-1.5,1.5]);
xlabel('x');
ylabel('y');
title('y=sin(10*\pi*x)/x');
[maxVal,maxIndex]=max(y);
plot(x(maxIndex),maxVal,'r*');
text(x(maxIndex),maxVal,{['x:' num2str(x(maxIndex))],['y:' num2str(maxVal)]});
[minVal,minIndex]=min(y);
plot(x(minIndex),minVal,'ro');
text(x(minIndex),minVal,{['x:' num2str(x(minIndex))],['y:' num2str(minVal)]});
hold off;

 

用模擬退火工具箱來找最值

求最小值

function fitness=fitnessfun(x)
fitness=sin(10*pi*x)./x;
end

 

 

 求最大值

function fitness=fitnessfun(x)
fitness=-sin(10*pi*x)./x;
end

Optimization running.
Objective function value: -0.9527670052175917
Maximum number of iterations exceeded: increase options.MaxIterations.
用工具箱求得的最大值為0.9527670052175917

 

 2.二元函數優化

[x,y]=meshgrid(-5:0.1:5,-5:0.1:5);
z=x.^2+y.^2-10*cos(2*pi*x)-10*cos(2*pi*y)+20;
figure
mesh(x,y,z);
hold on
xlabel('x');
ylabel('y');
zlabel('z');
title('z=x^2+y^2-10*cos(2*\pi*x)-10*cos(2*\pi*y)+20');
maxVal=max(z(:));
[maxIndexX,maxIndexY]=find(z==maxVal);%返回z==maxVal時，x和y的索引
for i=1:length(maxIndexX)
    plot3(x(maxIndexX(i),maxIndexY(i)),y(maxIndexX(i),maxIndexY(i)),maxVal,'r*');
    text(x(maxIndexX(i),maxIndexY(i)),y(maxIndexX(i),maxIndexY(i)),maxVal,{['x:' num2str(x(maxIndexX(i)))] ['y:' num2str(y(maxIndexY(i)))] ['z:' num2str(maxVal)] });
end
hold off;

 

function fitness=fitnessfun(x)
fitness=-(x(1).^2+x(2).^2-10*cos(2*pi*x(1))-10*cos(2*pi*x(2))+20);
end

 Optimization running.
Objective function value: -80.50038894455415
Maximum number of iterations exceeded: increase options.MaxIterations.
找到的最大值：80.50038894455415

 

 3.解TSP問題

（用的數據和前幾天用遺傳算法寫TSP問題的數據一致，但是結果比遺傳算法算出來效果差很多，不知道是不是我寫錯了，懷疑人生_(:з」∠)_中。。。

x=[41 94;37 84;54 67;25 62;7 64;2 99;68 58;71 44;54 62;83 69;64 60;18 54;22 60;
  83 46;91 38;25 38;24 42;58 69;71 71;74 78;87 76;18 40;13 40;82 7;62 32;58 35;45 21;41 26;44 35;4 50];
d=distance(x);%計算距離矩陣
n=size(d,1);%城市的個數
T0=1e50;%初始溫度
Tend=1e-30;%終止溫度
L=2;%各溫度下的迭代次數
q=0.95;
%計算迭代的次數 T0 * (0.95)^x = Tf
time=ceil(double(solve([num2str(T0) '*(0.9)^x=' num2str(Tend)])));%計算迭代的參數
count=0;
obj=zeros(time,1);
path=zeros(time,n);
%隨機產生一條初始路線
journey=randperm(n);
rlength=pathlength(d,journey);
drawpath(journey,x,rlength);
disp('初始種群中的一個隨機值:');
outputpath(journey);
disp(['初始路徑總距離：',num2str(rlength)]);
index=0;
%迭代優化
while T0>Tend
    count=count+1;
    %產生新解
    newjourney=NewAnswer(journey);
    [journey,r]=Metropolis(journey,newjourney,d,T0);
    %記錄每次迭代的最優路線
    if count==1||r<obj(count-1)
     obj(count)=r;
     index=count;
    else
        obj(count)=obj(count-1);
    end
     path(count,:)=journey;
     T0=T0*q;
     
     disp(['第' num2str(count) '代最優解' num2str(path(count,:))] );
     disp(['總距離：',num2str(obj(count))]);
end
    figure
    plot(1:count,obj);
    xlabel('迭代次數');
    ylabel('距離');
    title('優化過程');
grid on;
s=path(index,:);
a=pathlength(d,s);
drawpath(path(index,:),x,a);
disp('最優解');
p=outputpath(s);
disp(['總距離：',num2str(a)]);

 

function d=distance(citys)
%輸入 各城市的位置坐標
%輸出 兩兩城市之間的距離
n=size(citys,1);
d=zeros(n,n);
for i=1:n
    for j=1:n
        d(i,j)=((citys(i,1)-citys(j,1))^2+(citys(i,2)-citys(j,2))^2)^0.5;
        d(j,i)=d(i,j);
    end
end

 

function [s,r]=Metropolis(s1,s2,d,t)
%s1 當前解
%s2 新解
%d 距離矩陣
%t 溫度
%s 下一個當前解
%r 下一個當前解的路線距離
r1=pathlength(d,s1);%計算s1路線長度
n=size(d(2));%城市的個數
r2=pathlength(d,s2);%計算s2路線長度
dc=r2-r1;%計算能力之差

%比前一個解更優接受新解，沒有前一個解以一定概率接受新解
if dc<0 
    s=s2;
    r=r2;
elseif exp(-dc/t)>0 % 以exp(-dC/T)概率接受新路線
        a=exp(-dc/t)*100;
        test(1:100)=0;
        test(1:a)=1;
        r=round(99*rand+1);
        pcc=test(r);
        if pcc==1
        s=s2;
        r=r2;
        else
            s=s1;
            r=r1;
        end
    else
        s=s1;
        r=r1;
end
end

 

function s2=NewAnswer(s1)
%隨機產生新的解
%s1 當前解
%s2 新解
n=length(s1);
s2=s1;
a=round(rand(1,2)*(n-1)+1);
s2(a(2))=s1(a(1));
w=s1(a(2));
s2(a(1))=w;

 

function p=outputpath(route)
%給R末端添加起點R(1)
route=[route,route(1)];
n=length(route);
p=num2str(route(1));
for i=2:n
p=[p,'->',num2str(route(i))];
end
disp(p);

 

function dis=pathlength(d,route)
%計算起點與終點之間的距離
%d 城市間距離
%route 路線
dis=0;
n=length(d);
for i=1:n-1
    j=route(i);
    k=route(i+1);
   dis=dis+d(route(j),route(k));
   
end
j=route(n);
dis=dis+d(route(1),route(j));

 

function drawpath(route,citys,s)
%傳入參數 路線 城市坐標 
figure
plot([citys(route, 1); citys(route(1), 1)], [citys(route, 2); citys(route(1), 2)], 'o-');
grid on
 
%給每個地點標上序號
for i = 1: size(citys, 1)
    text(citys(i, 1), citys(i, 2), ['   ' num2str(i)]);
end
 
text(citys(route(1), 1), citys(route(1), 2), '       起點');
text(citys(route(end), 1), citys(route(end), 2), '       終點');
text(10,20,['距離長度：' num2str(s) 'km']);
xlabel('x/km');
ylabel('y/km');