概述

模拟退火的出发点是基于物理中固体物质的退火过程与一般的组合优化问题之间的相似性。

原理

1.加温过程

增强粒子的热运动，使其偏离平衡位置，当温度足够高的时候，固体将熔为液体，消除原有的非均匀状态。

2.等温过程

对于与周围环境交换热量而温度不变的封闭系统，系统状态的自发变化总是朝自由能减少的方向进行的，当自由能达到最小时，系统达到平衡状态。

3.冷却过程

使粒子热运动减弱，系统能量下降，得到晶体结构。

加温过程相当于对算法设定初值，等温过程对应算法的Metropolis抽样过程，冷却过程对应控制参数的下降。这里能量的变化就是目标函数，我们要得到的最优解就是能量最低态。其中Metropolis准则是SA算法收敛于全局最优解的关键所在，Metropolis准则以一定的概率接受恶化解，这样就使算法跳离局部最优的陷阱。

SA算法的Metropolis准则允许接受一定的恶化解，具体来讲，是以一定概率来接受非最优解。举个例子，相当于保留一些“潜力股”，使解空间里有更多的可能性。对比轮盘赌法，从概率论来讲，它是对非最优解给予概率0，即全部抛弃。

模拟退火本身是求一个最小值问题，但可以转化为求最大值问题，只需要对目标函数加个负号或者取倒数。

Metropolis准则

$P= \begin{cases} 1,&E_{t+1}<E_t\\ e^{\frac{-(E_{t+1}-E_t)}{kT}},&E_{t+1}\geq E_t \end{cases}$

特点

不属于群优化算法，不需要初始化。
收敛速度慢。
温度管理对结果影响大。

算法流程

v2-e1529b1b1577eecd16fb8a5418d51533_1440w

具体代码

非线性问题

clear 
clc
%%main.m
a=linspace(-10,10);
y=Fx(a);
T=1000; %初始化温度值
T_min=1e-12; %设置温度下界
alpha=0.98; %温度的下降率
k=1000; %迭代次数（解空间的大小）
 
x=getX; %随机得到初始解
while(T>T_min)
    for I=1:100
        fx=Fx(x);
        x_new=getX;
        if(x_new>=-2 && x_new<=2)
            fx_new=Fx(x_new);
            delta=fx_new-fx;
            if (delta<0)
                x=x_new+(2*rand-1);
            else
                P=getP(delta,T);
                if(P>rand)
                    x=x_new;
                end
            end
        end
    end
    T=T*alpha;
end
hold on
plot(a,y);
scatter(x,Fx(x))
disp('最优解为：')
disp(x)
hold off
%%getX.m
function x=getX
    x=4*rand-2;
end
 
%%Fx.m
function fx=Fx(x)
    fx=(x-2).^2+4;
end
 
%%getP.m
function p=getP(c,t)
    p=exp(-c/t);
end

旅行商问题

%% I. 清空环境变量
clear all
clc

%% II. 导入城市位置数据
X = [16.4700   96.1000
    16.4700   94.4400
    20.0900   92.5400
    22.3900   93.3700
    25.2300   97.2400
    22.0000   96.0500
    20.4700   97.0200
    17.2000   96.2900
    16.3000   97.3800
    14.0500   98.1200
    16.5300   97.3800
    21.5200   95.5900
    19.4100   97.1300
    20.0900   92.5500];

%% III. 计算距离矩阵
D = Distance(X);  %计算距离矩阵
n = size(D,1);    %城市的个数

%% IV. 初始化参数
T0 = 1e10;                                                                  % 初始温度
Tf = 1e-30;                                                                 % 终止温度
L = 2;                                                                      % 各温度下的迭代次数
q = 0.99;                                                                    % 降温速率
syms x;
eq1=T0*(0.99)^x ==Tf;
Time = ceil(double(solve(eq1,x)));                                          % 计算迭代的次数 T0 * (0.9)^x = Tf
count = 0;                                                                  % 初始化迭代计数
Obj = zeros(Time, 1);                                                       % 目标值矩阵初始化
path = zeros(Time, n);                                                      % 每代的最优路线矩阵初始化

%% V. 随机产生一个初始路线
S1 = randperm(n);
DrawPath(S1, X)                                                             % 画出初始路线
disp('初始种群中的一个随机值:')
OutputPath(S1);                                                             % 用箭头的形式表述初始路线
rlength = PathLength(D, S1);                                                % 计算初始路线的总长度
disp(['总距离：', num2str(rlength)]);

%% VI. 迭代优化
while T0 > Tf
    count = count + 1;                                                      % 更新迭代次数
    %temp = zeros(L,n+1);
    % 1. 产生新解
    S2 = NewAnswer(S1);
    % 2. Metropolis法则判断是否接受新解
    [S1, R] = Metropolis(S1, S2, D, T0);                                    % Metropolis 抽样算法
    % 3. 记录每次迭代过程的最优路线
    if count == 1 || R < Obj(count-1)                                       % 如果当前温度下的距离小于上个温度的距离，记录当前距离
        Obj(count) = R;
    else
        Obj(count) = Obj(count-1);
    end
    path(count,:) = S1;                                                     % 记录每次迭代的路线
    T0 = q * T0;                                                            % 以q的速率降温
end

%% VII. 优化过程迭代图
figure
plot(1: count, Obj)
xlabel('迭代次数')
ylabel('距离')
title('优化过程')
grid on

%% VIII. 绘制最优路径图
DrawPath(path(end, :), X)                                                   % path矩阵的最后一行一定是最优路线

%% IX. 输出最优解的路线和总距离
disp('最优解:')
S = path(end, :);
p = OutputPath(S);
disp(['总距离：', num2str(PathLength(D, S))]);
%% 计算两两城市之间的距离
function D = Distance(citys)

% 输入:各城市的位置坐标(citys)
% 输出:两两城市之间的距离(D)

n = size(citys, 1);
D = zeros(n, n);
for i = 1: n
    for j = i + 1: n
        D(i, j) = sqrt((citys(i, 1) - citys(j, 1))^2 + (citys(i, 2) - citys(j, 2))^2);
        D(j, i) = D(i, j);
    end
end
end

%% 画路径函数
function DrawPath(Route, citys)

% 输入
% Route  待画路径
% citys  各城市坐标位置

%画出初始路线
figure
plot([citys(Route, 1); citys(Route(1), 1)], [citys(Route, 2); citys(Route(1), 2)], 'o-');
grid on

%给每个地点标上序号
for i = 1: size(citys, 1)
    text(citys(i, 1), citys(i, 2), ['   ' num2str(i)]);
end

text(citys(Route(1), 1), citys(Route(1), 2), '       起点');
text(citys(Route(end), 1), citys(Route(end), 2), '       终点');
end
%% Metropolis准则判定
function [S,R] = Metropolis(S1, S2, D, T)

% 输入
% S1：  当前解
% S2:   新解
% D:    距离矩阵（两两城市的之间的距离）
% T:    当前温度
% 输出
% S：   下一个当前解
% R：   下一个当前解的路线距离

R1 = PathLength(D, S1);         % 计算S1路线长度
n = length(S1);                 % 城市的个数

R2 = PathLength(D,S2);          % 计算S2路线长度
dC = R2 - R1;                   % 计算能力之差
if dC < 0                       % 如果能力降低 接受新路线
    S = S2;                     % 接受新解
    R = R2;
elseif exp(-dC/T) >= rand       % 以exp(-dC/T)概率接受新路线
    S = S2;
    R = R2;
else                            % 不接受新路线
    S = S1;
    R = R1;
end
end
%% 产生新解
function S2 = NewAnswer(S1)
% 输入
% S1:    当前解
% 输出
% S2：   新解
 
n = length(S1);
S2 = S1;
a = round(rand(1, 2) * (n - 1) + 1);    %产生两个随机位置，用于交换
W = S2(a(1));
S2(a(1)) = S1(a(2));
S2(a(2)) = W;
end
%% 输出路径函数
function p = OutputPath(Route)

% 输入: 路径(R)
 
%给R末端添加起点即R(1)
Route = [Route, Route(1)];
n = length(Route);
p = num2str(Route(1));
 
for i = 2: n
    p = [p, ' —> ', num2str(Route(i))];
end
disp(p)
end
%% 计算起点与终点之间的路径长度
function Length = PathLength(D, Route)

% 输入：
% D     两两城市之间的距离
% Route 路线
 
Length = 0;
n = length(Route);
 
for i = 1: (n - 1)
    Length = Length + D(Route(i), Route(i + 1));
end
 
Length = Length + D(Route(n), Route(1));
end

参考文献

数学建模——模拟退火算法(Simulated Annealing，SA)_模拟退火法计算走时差-CSDN博客

模拟退火算法详解 - 知乎 (zhihu.com)

详解模拟退火算法（含MATLAB代码）_模拟退火 matlab-CSDN博客