模拟退火求解旅行商问题

Question

就是我代码是抄的，改了一点点，能够得出结果，但是结果跟我用lingo(也是抄的代码，lingo直接输个距离矩阵就能出结果)的相差甚远，并且这个模拟退火的路径也有问题，

clc,clear
%基于模拟退火求解
x=xlsread("附件1.xlsx",1,'B2:O15');

%模型假设
%可将无人机简化为
% 处理数据，将2点间无路线的距离设为10000；
DistMatrix = fillmissing(x,"constant",10000);
y_ = fillmissing(y,"constant",10000);%无人机的距离矩阵
%初始参数：T0(初始温度) T_end(结束温度) L(迭代次数) alpha(降温速率) S1(初始解) 
% 1.计算任意两点间的距离矩阵
% 2.生成新解
% 3.根据Metropolis准则判断是否接受新解
% 4.继续迭代，直到达到了迭代次数
% 5.降温
% 6.重复上面的2~5，直到温度降到最低
% 定义的函数：
% 1.根据每个城市的坐标计算任意距离矩阵函数Distence.m
% 2.绘制路线图Drawpath.m
% 3.计算一个解的距离的函数CalDist.m
% 4.生成新解的函数NewSolution.m
% 5.Metropolis准则函数,传入参数为(df,T_now),返回0或1。1代表接受新解,Metropolis.m
% 6.输出解的函数Disp.m
% 总体思路：给定初始温度，温度大于最低温度的时候，代表计算还在继续。每迭代一次，温度以指数规律减少。每次迭代时，指定链长L是操作的最大次数，我们每次迭代一共生成L次新解
%citydataset代表当前的(N*2)数据集
citydataset=x;
%模型的初始参数可以随时修改，L是链长，代表同一个温度下迭代的次数
T0=50000;T_end=1e-10;L=400;alpha=0.99;
[Num,~]=size(citydataset);%Num是城市的个数
S1=randperm(Num);%S1代表目前的解
disp(['初始路线的长度是：',num2str(CalDist(S1,DistMatrix))])
T_now=T0;
%Tset=T_now;
Tset=1;p=1;
Calset=CalDist(S1,DistMatrix);
while T_now>T_end
k=1;
while k<=L
S_new=NewSolution(S1);
df=CalDist(S_new,DistMatrix)-CalDist(S1,DistMatrix);
if Metropolis(df,T_now)==1 %接受新解
S1=S_new;
end
k=k+1;
end
T_now=T_now*alpha;
%Tset=[Tset,T_now];
Tset=[Tset,p];p=p+1;
cal=CalDist(S1,DistMatrix);
Calset=[Calset,cal];
end
Disp([S1])
disp(['最终路线的长度是：',num2str(CalDist(S1,DistMatrix))])
function L=Metropolis(df,T)
if df<0 %新解比当前解小的时候，一定接受新解
L=1;
else
r=rand(1);
if r<exp(-df/T) %新解大于当前解的时候，概率接受新解，接受概率取决于当前温度和新解与旧解之差
L=1;
else
L=0;
end
end
end
function sum=CalDist(s,Matrix)% s是当前解，Matrix是距离矩阵
[~,Num]=size(s);
sum=0;
for i=1:Num-1
sum=sum+Matrix(s(i),s(i+1));
end
sum=sum+Matrix(Num,1);
end
function s1=NewSolution(s)
[~,Num]=size(s);
R=round(rand(1,2)*(Num-1)+1);
s1=s;
s1(R(1))=s(R(2));
s1(R(2))=s(R(1));
end
function Disp(S)
[~,Num]=size(S);
p=num2str(S(1));
for i=2:Num-1
p=[p,'->',num2str(S(i))];
end
disp(p)
end>
 

下面是距离矩阵的图，感觉我代码没太大问题吧