数学建模——模拟退火算法(搜索全局最优解)

模拟退火算法

算法介绍

模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。

如果用粒子的能量定义材料的状态，Metropolis 算法用一个简单的数学模型描述了退火过程。假设材料在状态i 之下的能量为E(i)，那么材料在温度T 时从状态i 进入状态 j 就遵循如下规律：

1. 如果$E(j)\le E(i)$，则接受该状态被转化

2. 如果$E(j)>E(i)$，则该状态转换以如下概率被接受
   $$e^{\frac{E(i)-E(j)}{KT}}$$
   式中，K为物理学中的玻尔兹曼常数，T为材料温度

以下为算法的推导，不要求掌握，知道算法怎么用就可以

算法优点

模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。就是模拟退火算法能够有效避免陷入局部最优解，模拟退火能够帮我们有效的找到全局最优解。

例题

MATLAB代码详解

% 模拟退火算法
clc;
clear;
sj0=load('sj.txt');
x=sj0(:,[1:2:8]);x=x(:);
y=sj0(:,[2:2:8]);y=y(:);
% 观察数据存在的形式，转化成我们常见的形式
sj=[x,y];
dl=[70,40]; % 基地的位置
sj=[dl;sj;dl]; % 为了方便我们，我们将头和尾都设置为基地，中间为100个目标
sj=sj*pi/180;% 角度化成弧度
% 距离矩阵d初始化
d=zeros(102);
for i=1:101
    for j=i+1:102
        d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2)));
    end
end
% 根据上面推出的公式我们计算出每两点之间的距离，注意算出来的是矩阵对角线右上的数据，我们还得构造出实对称矩阵
d=d+d';% 实对称矩阵
path=[];long=inf;% 巡航路径及长度初始化,long为最远的距离，刚开始先初始化为无穷远
rand('state',sum(clock)); % 初始化随机模拟器
for j=1:1000 % 求较好的初始解，蒙特卡洛方法
    path0=[1 1+randperm(100),102];temp=0;
    % randperm是将[1:1:100]的顺序随机置换
    % path0是我们随机构造出的一个路径，头尾是基地，中间是100个目标的顺序
    for i = 1:101
        temp=temp+d(path0(i),path0(i+1));
        % 循环算出路径距离
    end
    if temp<long
        path=path0;
        long=temp;
        % 更新初始化的路径和距离
    end
end
e=0.1^30;L=20000;at=0.999;T=1;
% L退化次数
% T--初始温度
% at--每次退火的温度变化，即T'=T*0.999
% e--终止温度
for k=1:L % 退火过程
    c=2+floor(100*rand(1,2));% 产生新解
    c=sort(c);c1=c(1);c2=c(2);
    % 这里我们运用上面公式的2变化法，先随机选出两个u,v，即c(1)=u,c(2)=v
    % 计算代价函数增量
    df = d(path(c1-1),path(c2))+d(path(c1),path(c2+1))-d(path(c1-1),path(c1))-d(path(c2),path(c2+1));
    % 直接代上面的公式计算df
    if df<0 %接受准则
        path=[path(1:c1-1),path(c2:-1:c1),path(c2+1:102)];
        long=long+df;
    elseif exp(-df/T)>=rand
        path=[path(1:c1-1),path(c2:-1:c1),path(c2+1:102)];long=long+df;
    end
    T=T*at;
    if T<e
        break
    end
end
xx=sj(path,1);
yy=sj(path,2);
plot(xx,yy,'-*')% 画出巡航路径

巡航路径

附录

sj.txt

53.7121   15.3046	51.1758    0.0322	46.3253   28.2753	30.3313    6.9348
56.5432   21.4188	10.8198   16.2529	22.7891   23.1045	10.1584   12.4819
20.1050   15.4562	1.9451    0.2057	26.4951   22.1221	31.4847    8.9640
26.2418   18.1760	44.0356   13.5401	28.9836   25.9879	38.4722   20.1731
28.2694   29.0011	32.1910    5.8699	36.4863   29.7284	0.9718   28.1477
8.9586   24.6635	16.5618   23.6143	10.5597   15.1178	50.2111   10.2944
8.1519    9.5325	22.1075   18.5569	0.1215   18.8726	48.2077   16.8889
31.9499   17.6309	0.7732    0.4656	47.4134   23.7783	41.8671    3.5667
43.5474    3.9061	53.3524   26.7256	30.8165   13.4595	27.7133    5.0706
23.9222    7.6306	51.9612   22.8511	12.7938   15.7307	4.9568    8.3669
21.5051   24.0909	15.2548   27.2111	6.2070    5.1442	49.2430   16.7044
17.1168   20.0354	34.1688   22.7571	9.4402    3.9200	11.5812   14.5677
52.1181    0.4088	9.5559   11.4219	24.4509    6.5634	26.7213   28.5667
37.5848   16.8474	35.6619    9.9333	24.4654    3.1644	0.7775    6.9576
14.4703   13.6368	19.8660   15.1224	3.1616    4.2428	18.5245   14.3598
58.6849   27.1485	39.5168   16.9371	56.5089   13.7090	52.5211   15.7957
38.4300    8.4648	51.8181   23.0159	8.9983   23.6440	50.1156   23.7816
13.7909    1.9510	34.0574   23.3960	23.0624    8.4319	19.9857    5.7902
40.8801   14.2978	58.8289   14.5229	18.6635    6.7436	52.8423   27.2880
39.9494   29.5114	47.5099   24.0664	10.1121   27.2662	28.7812   27.6659
8.0831   27.6705	9.1556   14.1304	53.7989    0.2199	33.6490    0.3980
1.3496   16.8359	49.9816    6.0828	19.3635   17.6622	36.9545   23.0265
15.7320   19.5697	11.5118   17.3884	44.0398   16.2635	39.7139   28.4203
6.9909   23.1804	38.3392   19.9950	24.6543   19.6057	36.9980   24.3992
4.1591    3.1853	40.1400   20.3030	23.9876    9.4030	41.1084   27.7149