模拟退火算法matlab代码实现

采用模拟退火算法来实现路径优化。
即存在若干个不同位置的城市，如何安排旅程顺序，使从第一个城市旅行完所有城市并回到初始城市时所经过的距离最短。
首先放主程序demo.m

clc
clear
%%
tic
T0 = 1000;%初始温度
Tend = 1e-3;%终止温度
L = 200;%各温度下的迭代次数
q = 0.5;%降温速率
N= 20;%城市个数
X = rand(N,2)*10; %随机生成城市坐标
% X = [16.47 96.1
%      16.47 94.44
%      20.09 92.54
%      22.39 93.37
%      25.23 97.24
%      22 96.05
%      20.47 97.02
%      17.2 96.29
%      16.3 97.38
%      14.05 98.12
%      16.53 97.38
%      21.52 95.59
%      19.41 97.13
%      20.09 92.55];
D = Distanse(X);
S1 = randperm(N);%%随机生成初始解
S_initial = S1;%%随机生成初始解
Rlength = PathLength(D,S1);
disp(['初始值S：:',num2str(S1)])
disp(['初始值总距离：:',num2str(Rlength)])
syms x
% q = ['1000*(0.9)^x=',num2str(Tend)];
Time = ceil(double(solve(1000*q^x==Tend,x)));
count = 0;
Obj = zeros(Time,1);
track = zeros(Time,N);

while T0>Tend
    count = count+1;
    temp = zeros(L,N+1);
    for k = 1:L
        %%产生新解
        S2 = NewAnswer(S1);
        %%根据法则判断是否接受新解
        [S1,R] = Metropolis(S1,S2,D,T0);
        temp(k,:) = [S1,R];%记录每次迭代的旅游策略和对应的距离
    end
    %%记录每个温度下得到的最优解
    [d0,index] = min(temp(:,end));
    if count == 1 || d0<Obj(count-1)
        Obj(count) = d0;
    else
        Obj(count) = Obj(count-1);
    end

    track(count,:) = temp(index,1:end-1);%记录当前温度的最优路线
    T0 = q*T0;
    fprintf(1,'%d\n',count);
end
figure(1)
DrawPath(S_initial,X)%(X,S1)
title('初始值')
figure(2)
DrawPath(track(end,:),X)%(X,S1)
title('优化结果')
figure(3)
plot(1:Time,Obj)
title('优化过程')
%%输出结果
disp(['最优解：',num2str(track(end,:))]);
disp(['最优解对应的距离：',num2str(Obj(end))]);
toc

除此之外，还有几个函数的定义
函数1 NewAnswer.m用来生成新的策略

function S2 = NewAnswer(S1)
%NewAnswer 交换新解
% input： S1（上一步的解）
% output：S2（当前步的解）
N = length(S1);
S2 = S1;
a = randi(N,1,2);
if a(1)==a(2)
    a = randi(N,1,2);
end
W = S2(a(1));
S2(a(1)) = S2(a(2));
S2(a(2)) = W;
end

函数2 Distanse.m

function D = Distanse(a)
%Distanse 计算两两城市之间的距离
%input : a 
row = size(a,1);
D = zeros(row,row);
for i = 1:row
    for j = i:row
        D(i,j) = ((a(i,1)-a(j,1))^2+(a(i,2)-a(j,2))^2)^0.5;
        D(j,i) = D(i,j);
    end
end

函数3 PathLength.m

function len = PathLength(D,S)
%根据轨迹计算旅行总距离
%input ： D（记录城市两两之间距离的矩阵） Chorm 旅行城市的顺序
%output: len（旅行总距离） 

N = length(S);
len = 0;
for i = 1:N-1
    len = len + D(S(i),S(i+1));
    if i == N-1
        len = len+D(S(i+1),S(1));
    end
end

函数4 Metropolis.m 采用该方法判断是不是接受新解

function [S,R] = Metropolis(S1,S2,D,T)
%Metropolis 
%input : S1(当前解) S2(新解) D(距离矩阵) T(当前温度)
%output: S 下一个时刻的解 R(下一个时刻解的路线距离)
R1 = PathLength(D,S1);
N = length(S1);
R2 = PathLength(D,S2);
dC = R2-R1;
if dC<0
    S = S2;
    R = R2;
elseif exp(-dC/T)>=rand
    S = S2;
    R = R2;
else
    S = S1;
    R = R1;
end

函数5 DrawPath.m 绘制路线轨迹图（需要调用函数6）

function DrawPath(S,X)%(S,X)
%DrawPath 绘制路线图
%input : S(旅行策略) X(城市坐标)
% plot(1:10,1:10);
plot(X(:,1),X(:,2),'o','color',[0.5,0.5,0.5]);
for i = 1:size(X,1)
    text(X(i,1)+0.05,X(i,2)+0.05,num2str(i),'color',[1,0,0]);
end

for i =1:length(S)-1
    drawLine_Arrow(X(S(i),:),X(S(i+1),:),'r','g');  
    if i == length(S)-1
        drawLine_Arrow(X(S(i+1),:),X(S(1),:),'r','g');
    end
end
end

函数6 drawLine_Arrow.m

function drawLine_Arrow(start_point, end_point,arrColor,lineColor,arrowSize,lineWidth)
%% 绘制带箭头的直线
% 从start_point到end_point画一箭头,arrColor箭头颜色,arrSize，箭头大小
%判断参数多少
switch nargin
    case 2
        arrColor  = 'r';
        lineColor = 'b';
        arrowSize = 12;
        lineWidth = 1;
    case 3
        lineColor = 'b';
        arrowSize = 12;
        lineWidth = 1;
    case 4
        arrowSize = 12;
        lineWidth = 1;
    case 5
        lineWidth = 1;
end
K= 0.05;          % 箭头比例系数
theta= pi / 8;    % 箭头角度
% theta         =atan2(end_point(1,2)-start_point(1,2),end_point(1,1)-start_point(1,1));
% arrow_1(1)    = end_point(1,1)-arrowSize*cos(theta)+K*sin(theta);
% arrow_2(1)    = end_point(1,1)-arrowSize*cos(theta)-K*sin(theta);
% arrow_1(2)    = end_point(1,2)-arrowSize*sin(theta)-K*cos(theta);
% arrow_2(2)    = end_point(1,2)-arrowSize*sin(theta)+K*cos(theta);

A1 = [cos(theta), -sin(theta);sin(theta), cos(theta)]; % 箭头左侧
theta = -theta;
A2 = [cos(theta), -sin(theta);sin(theta), cos(theta)]; % 箭头右侧
arrow = start_point' - end_point'
%使得箭头跟直线长短无关(固定值)
% arrow(arrow>=0) = -arrowSize 
% arrow(arrow<0) = -arrowSize
arrow_1(1)    = end_point(1,1)-arrowSize*cos(theta)+K*sin(theta);
arrow_2(1)    = end_point(1,1)-arrowSize*cos(theta)-K*sin(theta);
arrow_1(2)    = end_point(1,2)-arrowSize*sin(theta)-K*cos(theta);
arrow_2(2)    = end_point(1,2)-arrowSize*sin(theta)+K*cos(theta);

arrow_1= A1 * arrow;
arrow_2= A2 * arrow;
arrow_1= K * arrow_1 + end_point'; % 箭头的边的x坐标
arrow_2= K * arrow_2 + end_point'; % 箭头的边的y坐标
hold on;
grid on;
axis equal;
plot([start_point(1), end_point(1)], [start_point(2), end_point(2)],lineColor,'lineWidth',lineWidth);
% 三角箭头(填充)
triangle_x= [end_point(1),arrow_1(1),arrow_2(1),end_point(1)];
triangle_y= [end_point(2),arrow_1(2),arrow_2(2),end_point(2)];
fill(triangle_x,triangle_y,arrColor);
% 线段箭头(不填充)
% plot([arrow_1(1), end_point(1)], [arrow_1(2), end_point(2)],color,'lineWidth',arrowSize);
% plot([arrow_2(1), end_point(1)], [arrow_2(2), end_point(2)], color,'lineWidth',arrowSize);
hold on;
end

结果展示
1、初始轨迹

2、最终优化之后的轨迹

可以看出优化后的轨迹更加合理

3、距离随迭代次数收敛情况

横轴是迭代次数，纵轴是旅行总距离，这里忘了标了。

参考资料：
1、MATLAB智能算法30个案例分析
2、Matlab工具— 为两点绘制带箭头的连线（函数6来源）