load ij.txt;%加載目标資料ij.txt存放矩陣
x=ij(:,1:2:8);
%把sj中屬于1,3,5,7的列所有元素賦給x,此時,x的行的大小為sj的行的大小,列的大小為4,
x=x(:);%x=x(:)的意思是把x拉成一維列向量(應該是列向量)
y=ij(:,2:2:8);
y=y(:);
sj=[x y];%将xy合并成一個矩陣
d1=[70,40];%起始點
sj=[d1;sj;d1];%首尾添加上出發點和終點
sj=sj*pi/180;%角度化為弧度
%計算距離
d=zeros(102);%先建一個102x102的空的距離矩陣，再把對應距離填進去
%A(R cos x1 cos y1 , R sin x1 cos y1 , R sin y1 )
%B(R cos x2 cos y2 , R sin x2 cos y2 , R sin y2 )
%其中 R = 6370 為地球半徑。
%A, B 兩點的實際距離:
%d = R arccos[cos(x1 − x2 ) cos y1 cos y2 + sin y1 sin y2 ]
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=d+d';
S0=[];
Sum=inf;
rand('state',sum(clock));%作用是定義一個随時間變化的初值
%matlab裡面的随機生成函數基本都是以rand為基函數通過函數關系式得到，
%比如 normrnd，unidrnd等，你每次重新開機matlab後運作已編好的含随機數生成的函數你将得到相同的結果，
%比如我的電腦上重新開機matlab運作unidrnd(100),每次的值都是82，這是因為rand函數的初值都一樣，
%是以為了避免上述問題經常在程式前運作或加指令rand('state',sum(clock))，這樣重新開機matlab，運作随機數生成值就不同了。
for i=1:1000%循環1000次，盡可能讓所有排序情況都出現一遍
s=[1,1+randperm(100),102];%s為中間2到101的排列組合
temp=0;
for j=1:101
temp=temp+d(s(j),s(j+1));
end
if temp<Sum
Sum=temp;
S0=s;%儲存目前最小值的排序情況
end
end
e=0.1^30;%標明的終止溫度 e = 10 −30 ，判斷退火過程是否結束
at=0.999;%標明的降溫系數α進行降溫即：T =αT
L=20000;
T=1;
%退火過程
for k=1:L
c=2+floor(100*rand(1,2));%floor函數其功能是“向下取整”，或者說“向下舍入”，即取不大于x的最大整數
%ceil函數向上取整
c=sort(c);% [a,b]=sort(X)是按列從小到大排序，而 [a,b]=sort(X,2)是按行，b為排序情況
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));
if df<0%接受準則
S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
Sum=Sum+df;
elseif exp(-df/T)>rand(1)%以機率exp(−df/T)接受新的路徑%注意時elseif而不是else if
S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
Sum=Sum+df;
end;
T=T*at;%退火
if T<e%達到終止溫度
break;
end;
end;
%輸出路徑和路徑長度
S0,Sum
%畫出其中一個巡航路徑
plot(sj(S0,1)/pi*180,sj(S0,2)/pi*180);
hold on
plot(sj(S0,1)/pi*180,sj(S0,2)/pi*180,'rx');
axis([-5,75,-5,45]);
![數學模組化淺談[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)