Posted by on 2017年3月14日

2 comments

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

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

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

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

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

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

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

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

MATLAB
1
2
3
>> len(1)
ans =
   1.0855e+03

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

MATLAB
1
2
3
4
5
6
7
8
9
10
11
% 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如下：

MATLAB
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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
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

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

MATLAB
1
2
3
4
5
6
7
8
9
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

画图函数如下：

MATLAB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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

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

MATLAB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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

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

总路长为：

MATLAB
1
2
3
>>len(end)
ans =
  406.3937

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

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

Posted in: Academic 無涯齋
Tags: , ,
120 views
1 Star2 Stars3 Stars4 Stars5 Stars (No Ratings Yet)
Loading...

Comments

  1. luo
    2017年3月15日

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

Leave a Reply

返回顶部