一个模拟退火(SA)的小实验

Posted by baoduge on 2017年3月14日

今天来做个模拟退火(Simulated Annealing，百度百科点我！，wikipedia click here!)的小实验。

模拟退火是一种通用概率算法，用来在固定时间内寻求在一个大的搜寻空间内找到的最优解。模拟退火是S. Kirkpatrick, C. D. Gelatt和M. P. Vecchi在1983年所发明。而V. Černý在1985年也独立发明此算法。

我觉得它好玩的地方就在于该算法模拟了冶金中的退火过程，同时也让我重新认识到材料的退火其实也是一个最优化的过程。退火通过反复地升温和降温（淬火），不断地寻找分子化学式的全局最优解，同时也避免最终状态局限在局部最优。

模拟退火算法与初始值无关，算法求得的解与初始解状态S（是算法迭代的起点）无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率1收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。

下面我们来做个简单的小实验吧。

我们的问题是：我们要规划一个最短的旅行路线可以达到地图上的所有城市。

首先我们随机生成一个地图如下：

地图上有20个城市，我们对他们随机进行标号，假设我们的初始路线是按照标号从1走到n=20，并计算初始总路长：

>> len(1)
ans =
   1.0855e+03

>> len(1)

ans =

1.0855e+03

下面我们设置模拟退火的参数如下：

% We are planning a trip with the shortest length and reaching eavry city
% in the map
n=20;                   % # of cities
temperature=100*n;      % initial temperature
iter=100;              
% evaluate the total # of iteration of every time of annealing
precision=1e-3;         % tolerance of termination
Vol_Decline=0.99;       % Rate of cooling down

% We are planning a trip with the shortest length and reaching eavry city

% in the map

n=20; % # of cities

temperature=100*n; % initial temperature

iter=100;

% evaluate the total # of iteration of every time of annealing

precision=1e-3; % tolerance of termination

Vol_Decline=0.99; % Rate of cooling down

接下来就是最主要的模拟退火过程：

1.随着温度不断下降，我们每次随机调换两个城市的标号，并且计算前后的总路长。

2.如果路长变短则采用调整后的规划路线，但是有一定的概率不采取更短的路线，而是保持原有的更长的路线（“淬火”）。概率模型为 exp(-两次前后路程差/当时温度)，温度越低，“淬火”概率越大；

3.每次调整后，温度都会以一定比率降低;

4.当温度达到容忍值时，终止程序.

具体的code如下：

N_Cal=log(precision/temperature)/log(Vol_Decline);
% Total # of annealing

% generate random map of cities
% suppose that the x,y coordinate is within [0,100]
city=struct([]);
% the initial plan is travel from 1 to n
for i=1:n
    city(i).x=100*rand(); 
    city(i).y=100*rand();
end

l=1;                            % # of annealing
len(l)=computer_tour(city,n);   % total length of trip
netplot(city,n);                % plot the initial planning
saveas(gcf,'Initial Tour','png');
savefig(gcf,'Initial Tour');

bar=waitbar(0,'Please Wait...:0%');% start a waiting bar
while temperature>precision    % terminating condition
    percentage=floor(1+100*l/N_Cal)/100;% percentage of computation process
    waitbar(percentage,bar,strcat('Please Wait:',num2str(100*percentage),'%'));
    for i=1:iter     
        len1=computer_tour(city,n);         
        tmp_city=perturb_tour(city,n);      % reorder my trip randomly
        len2=computer_tour(tmp_city,n);     % 
        
        delta_e=len2-len1;  
        % difference of the length in continuous computation
        if delta_e<0        % total length reduce
            city=tmp_city;
        else
            % there is possibility the total length will grow
            if exp(-delta_e/temperature)>rand() % 
                city=tmp_city;      % 
            end
        end        
    end
    l=l+1;
    len(l)=computer_tour(city,n);   % compute the total length of trip
    temperature=temperature*Vol_Decline;   
    % after every time of annealing 
  
end 
close(bar);
beep

figure;
netplot(city,n);    % plot my final planning of trip
saveas(gcf,'Optimized Tour','png');
savefig(gcf,'Optimized Tour');
figure;
plot(len);xlabel('Iteration');ylabel('Length');title('Length of City Tour');
saveas(gcf,'Length of Tour','png');
savefig(gcf,'Length of Tour');
toc

N_Cal=log(precision/temperature)/log(Vol_Decline);

% Total # of annealing

% generate random map of cities

% suppose that the x,y coordinate is within [0,100]

city=struct([]);

% the initial plan is travel from 1 to n

for i=1:n

city(i).x=100*rand();

city(i).y=100*rand();

end

l=1; % # of annealing

len(l)=computer_tour(city,n); % total length of trip

netplot(city,n); % plot the initial planning

saveas(gcf,'Initial Tour','png');

savefig(gcf,'Initial Tour');

bar=waitbar(0,'Please Wait...:0%');% start a waiting bar

while temperature>precision % terminating condition

percentage=floor(1+100*l/N_Cal)/100;% percentage of computation process

waitbar(percentage,bar,strcat('Please Wait:',num2str(100*percentage),'%'));

for i=1:iter

len1=computer_tour(city,n);

tmp_city=perturb_tour(city,n); % reorder my trip randomly

len2=computer_tour(tmp_city,n); %

delta_e=len2-len1;

% difference of the length in continuous computation

if delta_e<0 % total length reduce

city=tmp_city;

else

% there is possibility the total length will grow

if exp(-delta_e/temperature)>rand() %

city=tmp_city; %

end

l=l+1;

len(l)=computer_tour(city,n); % compute the total length of trip

temperature=temperature*Vol_Decline;

% after every time of annealing

end

close(bar);

beep

figure;

netplot(city,n); % plot my final planning of trip

saveas(gcf,'Optimized Tour','png');

savefig(gcf,'Optimized Tour');

figure;

plot(len);xlabel('Iteration');ylabel('Length');title('Length of City Tour');

saveas(gcf,'Length of Tour','png');

savefig(gcf,'Length of Tour');

toc

其中
计算总路长的函数如下：

function len=computer_tour(city,n)   
% computing the total length of trip
len=0;
for i=1:n-1
    len=len+sqrt((city(i).x-city(i+1).x)^2+(city(i).y-city(i+1).y)^2);        
end
len=len+sqrt((city(n).x-city(1).x)^2+(city(n).y-city(1).y)^2);
end

function len=computer_tour(city,n)

% computing the total length of trip

len=0;

for i=1:n-1

len=len+sqrt((city(i).x-city(i+1).x)^2+(city(i).y-city(i+1).y)^2);

end

len=len+sqrt((city(n).x-city(1).x)^2+(city(n).y-city(1).y)^2);

end

画图函数如下：

function netplot(city,n)        
% plot the map
    hold on;
    for i=1:n-1
        plot(city(i).x,city(i).y,'r*'); 
        text(city(i).x+0.1,city(i).y+0.1,strcat('City:',num2str(i)));
        line([city(i).x city(i+1).x],[city(i).y city(i+1).y]);  %????????????     
    end

    plot(city(n).x,city(n).y,'r*'); 
    text(city(n).x+0.1,city(n).y+0.1,strcat('City:',num2str(n)));
    line([city(n).x city(1).x],[city(n).y city(1).y]);     %?????????????
    hold off;
    xlabel('x');ylabel('y');title('Tour Plot(SA)');
end

function netplot(city,n)

% plot the map

hold on;

for i=1:n-1

plot(city(i).x,city(i).y,'r*');

text(city(i).x+0.1,city(i).y+0.1,strcat('City:',num2str(i)));

line([city(i).x city(i+1).x],[city(i).y city(i+1).y]); %????????????

end

plot(city(n).x,city(n).y,'r*');

text(city(n).x+0.1,city(n).y+0.1,strcat('City:',num2str(n)));

line([city(n).x city(1).x],[city(n).y city(1).y]); %?????????????

hold off;

xlabel('x');ylabel('y');title('Tour Plot(SA)');

end

随机城市排列的函数如下：

function city=perturb_tour(city,n)  
    p1=floor(1+n*rand());
    p2=floor(1+n*rand());
    % randomly pick up two cities
    while p1==p2 
        % if I picked two same cities
        p1=floor(1+n*rand());
        p2=floor(1+n*rand());    
    end
    
    tmp=city(p1);
    % exchange the order of this two 
    city(p1)=city(p2);
    city(p2)=tmp;
end

function city=perturb_tour(city,n)

p1=floor(1+n*rand());

p2=floor(1+n*rand());

% randomly pick up two cities

while p1==p2

% if I picked two same cities

p1=floor(1+n*rand());

p2=floor(1+n*rand());

end

tmp=city(p1);

% exchange the order of this two

city(p1)=city(p2);

city(p2)=tmp;

end

最终我们的规划路线如下图

总路长为：

>>len(end)
ans =
  406.3937

>>len(end)

ans =

406.3937

其中模拟退火计算过程中的总路长动态如下：

可以看出模拟退火（SA）算法的缺点也很明显，收敛速度慢，执行时间长。当然，计算最短路径的图算法还有很多，诸如Floyd-Warshall、Dijkstra‘s、Bellman-Ford、SPFA(Shortest Path Faster Algorithm）等，下次有机会再写。
【欢迎订阅并关注微信公众号“饱蠹阁BaoDuGe”】

Posted in: Academic 無涯齋

Tags: Algorithm, MATLAB, ToyModel

82 views

(No Ratings Yet)

Comments

luo
2017年3月15日

登录以回复

觉得很实用，举例也很清晰，利用退火思想在一点时间内，求一个全局最优解而不是局部最优解，就像退火的目的就是为了使铸件内分子结构更加均匀化，而达到提高整个铸件的物理性能。很棒！！！谢谢楼主讲解。
- baoduge
  2017年3月18日
  
  登录以回复
  
  谢谢评论

一	二	三	四	五	六	日
« 12月
1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28
29	30

一个模拟退火(SA)的小实验

Comments

Leave a Reply Cancel reply

微信公众号

Search

Recent Post

功能

一个模拟退火(SA)的小实验

Comments

Leave a Reply Cancel reply

微信公众号

Search

Labels

Recent Post

功能