模拟退火例题

例已知敌方100个目标的经度、纬度如下：经度纬度经度纬度经度纬度经度纬度53.712115.304651.17580.032246.325328.275330.33136.934856.543221.418810.819816.252922.789123.104510.158412.481920.105015.45621.94510.205726.495122.122131.48478.964026.241818.176044.035613.540128.983625.987938.472220.173128.269429.001132.19105.869936.486329.72840.971828.14778.958624.663516.561823.614310.559715.117850.211110.29448.15199.532522.107518.55690.121518.872648.207716.888931.949917.63090.77320.465647.413423.778341.86713.566743.54743.906153.352426.725630.816513.459527.71335.070623.92227.630651.961222.851112.793815.73074.95688.366921.505124.090915.254827.21116.20705.144249.243016.704417.116820.035434.168822.75719.44023.920011.581214.567752.11810.40889.555911.421924.45096.563426.721328.566737.584816.847435.66199.933324.46543.16440.77756.957614.470313.636819.866015.12243.16164.242818.524514.359858.684927.148539.516816.937156.508913.709052.521115.795738.43008.464851.818123.01598.998323.644050.115623.781613.79091.951034.057423.396023.06248.431919.98575.790240.880114.297858.828914.522918.66356.743652.842327.288039.949429.511447.509924.066410.112127.266228.781227.66598.083127.67059.155614.130453.79890.219933.64900.39801.349616.835949.98166.082819.363517.662236.954523.026515.732019.569711.511817.388444.039816.263539.713928.42036.990923.180438.339219.995024.654319.605736.998024.39924.15913.185340.140020.303023.98769.403041.108427.7149我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为1000公里/小时。我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。这是一个旅行商问题。我们依次给基地编号为1，敌方目标依次编号为2，3，…，101，最后我方基地再重复编号为102（这样便于程序中计算）。距离矩阵102102)(ijdD，其中ijd表示表示ji,两点的距离，102,,2,1,ji，这里D为实对称矩阵。则问题是求一个从点1出发，走遍所有中间点，到达点102的一个最短路径。上面问题中给定的是地理坐标（经度和纬度），我们必须求两点间的实际距离。设BA,两点的地理坐标分别为),(11yx，),(22yx，过BA,两点的大圆的劣弧长即为两点的实际距离。以地心为坐标原点O，以赤道平面为XOY平面，以0度经线圈所在的平面为XOZ平面建立三维直角坐标系。则BA,两点的直角坐标分别为：)sin,cossin,coscos(11111yRyxRyxRA)sin,cossin,coscos(22222yRyxRyxRB其中6370R为地球半径。BA,两点的实际距离OBOBRdOAOAarccos，化简得]sinsincoscos)(arccos[cos212121yyyyxxRd。求解的模拟退火算法描述如下：（1）解空间解空间S可表为{102,101,,2,1}的所有固定起点和终点的循环排列集合，即}102,}101,,3,2{),,(,1|),,{(102101211021的循环排列为S其中每一个循环排列表示侦察100个目标的一个回路，ji表示在第i次侦察j点，初始解可选为)102,,2,1(，本文中我们使用MonteCarlo方法求得一个较好的初始解。（2）目标函数此时的目标函数为侦察所有目标的路径长度或称代价函数。我们要求1011211102),,,(miniiidf而一次迭代由下列三步构成：（3）新解的产生①2变换法任选序号vu,（vu）交换u与v之间的顺序，此时的新路径为：10211111vuuvvu②3变换法任选序号vu,和w，将u和v之间的路径插到w之后，对应的新路径为（设wvu）1021111wvuwvu（4）代价函数差对于2变换法，路径差可表示为)()(1111vvuuvuvuddddf（5）接受准则0)/exp(01fTffP如果0f，则接受新的路径。否则，以概率)/exp(Tf接受新的路径，即若)/exp(Tf大于0到1之间的随机数则接受。（6）降温利用选定的降温系数进行降温即：TT，得到新的温度，这里我们取9999.0。（7）结束条件用选定的终止温度3010e，判断退火过程是否结束。若eT，算法结束，输出当前状态。我们编写如下的matlab程序如下：clc,clearloadsj.txt%加载敌方100个目标的数据，数据按照表格中的位置保存在纯文本文件sj.txt中x=sj(:,1:2:8);x=x(:);y=sj(:,2:2:8);y=y(:);sj=[xy];d1=[70,40];sj=[d1;sj;d1];sj=sj*pi/180;%距离矩阵dd=zeros(102);fori=1:101forj=i+1:102temp=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);endendd=d+d';S0=[];Sum=inf;rand('state',sum(clock));forj=1:1000S=[11+randperm(100),102];temp=0;fori=1:101temp=temp+d(S(i),S(i+1));endiftempSumS0=S;Sum=temp;endende=0.1^30;L=20000;at=0.999;T=1;%退火过程fork=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));%接受准则ifdf0S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];Sum=Sum+df;elseifexp(-df/T)rand(1)S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];Sum=Sum+df;endT=T*at;ifTebreak;endend%输出巡航路径及路径长度S0,Sum计算结果为44小时左右。其中的一个巡航路径如下图所示：0102030405060700510152025303540