智能算法之模拟退火算法

• • 1.起源
  • 2.物理退火流程
  • • 2.1 加温过程
    • 2.2 等温过程
    • 2.3 冷却过程
    • 2.4 组合优化与物理退化
  • 3.原理
  • • 3.1 算法核心迭代
    • 3.2 具体流程
  • 4.案例
  • • 4.1 求解n元函数的极小值
    • 4.2 求解二元函数的极小值
    • 4.3 货物配送规划
    • • 4.3.1 分析
    • 4.4 TSP问题求解
    • 4.5 车辆路径优化
  • 5.优缺点及可改进方向
  • • 5.1 优点
    • 5.2 缺点
    • 5.3 可改进方向
  • 6.参考文献

1.起源

模拟退火算法来源于热力学中固体物质的退火冷却过程 (退火是一种金属热处理工艺，指的是将金属缓慢加热到一定温度，保持足够时间，然后以适宜速度冷却)。

2.物理退火流程

模拟退火的算法思想是参考物理退火的过程而来，物理退火的过程为：加温过程 -> 等温过程 -> 冷却过程。

2.1 加温过程

通过加温，增强粒子的热运动，让各部分变得均匀。当温度足够高时，固体将熔解为液体，从而消除原先系统中可能存在的非均匀状态，保证后续进行的过程是从平衡态开始，由于是加温过程温度上升，所以该过程系统能量增加。

2.2 等温过程

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

2.3 冷却过程

通过一定速率对液体进行冷却，减弱粒子的热运动并趋于有序，在每个温度都达到平衡态，最后在常温时达到基态，能量减为最小，从而得到低能态的晶体结构。

2.4 组合优化与物理退化

组合优化	物理退火
解	粒子状态
最优解	能量最低态
设定初温	熔解过程
$Metropolis$ 抽样过程	等温过程
控制参数下降	冷却过程
目标函数	能量

3.原理

模拟退火算法的基本原理：对于目标函数 $f(X)$、自变量为 $X$ 的 $n$ 维极小化问题（极大化问题乘以$-1$转化为极小化问题），初值 $X_0$ 可以随机化或者自己设置。

3.1 算法核心迭代

设 $f_k,f_{k+1}$ 分别为目标函数在第 $k$ 次和第 $k+1$ 次迭代值，即 $f_k=f(X_k),f_{k+1}=f(X_{k+1})$。
若 $f_k>f_{k+1}$，则接受 $X_{k+1}$ 作为下一次迭代的初值继续进行运算，直至满足终止条件（收敛或者迭代次数耗尽）；
若 $f_k<f_{k+1}$，则接受 $X_{k+1}$ 的概率为 $p$，拒绝的概率为 $1-p$，其中接受概率公式如下：$$p=exp(-\frac{f_{k+1}-f_k}{T})$$ 其中 $T$ 为控制参数，在模拟退火算法中，$T$ 必须缓慢减少，避免变化过快，导致优化陷入极值点。

3.2 具体流程

以求解极小值问题为例：
$（1）$设定初始温度 $T$，迭代次数 $L$，误差 $ϵ$，初始解 $X_0$，得到目标函数$f(X_0)$；
$（2）$循环进行迭代，$k=1,\cdots ,L$；
$（3）$根据当前温度得到新解 $X^{\prime }$(第一个可以随机生成)， 代入 $f(X)$ 得到 $f(X^{\prime })$；
$（4）$若 $f(X)>f(X^{\prime })$，则接受 $X^{\prime }$ 作为下一次迭代的初值，若 $f(X)<f(X^{\prime })$，则以概率 $exp(-\frac{f(X^{\prime })-f(X)}{T})$ 接受 $X^{\prime }$ 作为下一次迭代的初值；
$（5）$若满足终止条件，则输出当前最优解 $X^{\prime }$，终止条件；
$（6）$$T$ 减小，然后回到流程$（2）$，在当前温度下继续循环迭代，执行步骤$（3，4，5）$；
$（7）$注意，这里的终止条件可以有多个，比如：迭代次数耗尽、满足最优解、温度降低到一定值、新目标函数与原先目标函数的误差小于 $ϵ$。

4.案例

4.1 求解n元函数的极小值

例1：求解函数 $f(X)=X^2$ 的最小值，其中 $X^2={\sum }_{i=1}^{10}{x}^2$，$-10\le x\le 10$。
下面通过matlab进行求解：
目标函数：

%--------目标函数f(x)---------
function result = func1(x)
    result = sum(x.^2);
end

模拟退火算法：

clear; clc;

%---------------------模拟退火核心-------------------------------
L = 200;                           % 每个温度的迭代次数
epsilon = 1e-10;                   % 新目标函数与原先目标函数的误差
T0 = 200;                          % 初始温度
T1 = 0.00001;                      % 终止温度
s = 0.01;                          % 自变量更新的步长，相当于学习率
K = 0.998;                         % 温度衰减参数
P = 0;                             % Metropolis 过程中总接受点

%---------------------根据题目所写--------------------------------
opt_minmax = 1;                    % 求解极小值问题为1，极大值问题为-1
n = 10;                            % 自变量维度
up = 10;                           % 自变量上限
sup = -10;                         % 自变量下限
x0 = rand(n, 1)*(up-sup)-up;       % 随机生成在范围内的初始自变量
x1 = rand(n, 1)*(up-sup)-up;       % 随机生成在范围内的 x1 自变量进行第一步
PreBest = x0;
NowBest = x1; 

%---------------------进行迭代----------------------------------
delta = abs(func1(NowBest) - func1(PreBest));
while (delta > epsilon) && (T0 > T1)
    T0 = K*T0;
    for i = 1:L
        x2 = x1 + s*(rand(n, 1)*(up-sup)-up);
        %-----下面这个循环是写题目的条件的，这里没啥好些，就重复了边界条件
        for j = 1:n
            while (x2(j) > up) || (x2(j) < sup)
                x2(j) = x1(j) + s*(rand*(up-sup)-up);
            end
        end
        %--------------查看解是否是全局最优解----------------
        if opt_minmax*func1(NowBest) > opt_minmax*func1(x2)
            PreBest = NowBest;             % 保留上一个最优解
            NowBest = x2;                  % 保留新的最优解
        end
        %---------------- Metropolis过程--------------------
        if opt_minmax*func1(x1) > opt_minmax*func1(x2)
            x1 = x2;
            P = P + 1;
        else
            p = exp(-(func1(x2)-func1(x1))/T0);
            if p > rand(1)
                x1 = x2;
                P = P + 1;
            end
        end
        Y(P) = func1(NowBest);
    end
    delta = abs(func1(NowBest) - func1(PreBest));
end

得到结果为：
$X=[-0.0114,-0.0002,0.0020,0.0101,0.0074,-0.0003,-0.0084,0.0112,-0.0252,0.0052]$

$f(X)=0.0011$
绘制适应度进化曲线：

%----------------------------画图--------------------------------
plot(Y, 'r');
hold on
scatter(length(Y), Y(end), 'bo', 'MarkerFaceColor', 'b')
text(length(Y), Y(end), ' 最优点')
xlabel('迭代次数'); ylabel('目标函数值'); title('适应度进化曲线');

4.2 求解二元函数的极小值

例2：求解函数 $f(x,y)=5cos(xy)+xy+y^3$ 的最小值，其中 $-5\le x\le 5,-5\le y\le 5$。
下面通过matlab进行求解：
目标函数：

%--------目标函数f(x)---------------
function result = func1(x)
    result = 5*cos(x(1).*x(2))+x(1).*x(2)+x(2).^3;
end

可以看看 $f(x,y)$ 的图像：

clear;clc;
x = linspace(-5, 5);
y = linspace(-5, 5);
[x, y] = meshgrid(x, y);
f = 5*cos(x.*y)+x.*y+y.^3;
mesh(x, y, f);
xlabel('x'); ylabel('y'); zlabel('f');

第二种方法适合小白，不过和meshgrid画出来的图相比较发现是翻转过的，所以可以自行选择。

x1 = linspace(-5, 5);
y1 = linspace(-5, 5);
N = size(x1, 2);
for i = 1:N
    for j = 1:N
        z(i, j)=5*cos(x1(i)*y1(j))+x1(i)*y1(j)+y1(j)^3;
    end
end
mesh(x1, y1, z);
xlabel('x1'); ylabel('y1'); zlabel('z');

这次使用matlab自带的模拟退火函数进行求解：

clear; clc;
%---------------------根据题目所写---------------------------------------
options = optimoptions('simulannealbnd','PlotFcns',...
          {@saplotbestx,@saplotbestf,@saplotx,@saplotf});
up = [5, 5];                         % 自变量x,y上限
sup = [-5, -5];                      % 自变量x,y下限
x0 = rand(1, 2).*(up-sup)-up;        % 随机生成在范围内的初始自变量x0,y0
func = @(x)(5*cos(x(1)*x(2))+x(1)*x(2)+x(2)^3);
[x, fval] = simulannealbnd(func, x0, sup, up, options)

得到结果为： $$f(x,y)=f(4.4242,-5.0000)=-152.0788$$

4.3 货物配送规划

例3：根据下面的所给的条件，给此地的生鲜配送路线进行规划，让生鲜被高效率送到各销售市场中，使得资源达到最大的利用。

$n$ ：销售市场数量；
$q_i$：每个市场需求量 $(i=1,2,\cdots ,n)$；
$m$：配送车数量（每辆车都型号都一致）；
$Q$：配送车最大装载量；
$d_{ij}$：市场 $i$ 到市场 $j$ 的距离；
$t_{ij}$：从市场 $i$ 到市场 $j$ 所需时间；
$d_{oi}$：配送中心 $o$ 到市场 $i$ 的距离；
$c_0$ ：每次配送固定成本；
$c_v$：单位里程行驶费用；
$v$：配送车匀速行驶速度；
$T$：生鲜装卸时间；
${\beta }_1$：单位时间配送运输损耗比例；
${\beta }_2$：单位时间配送装卸损耗比例；
$u_{ijk}$：配送车辆 $k$ 在 $i$ 到 $j$ 路径上行驶时的载重量；
$p$：单位重量货物的货损价格；
$f_0$：单次配送任务中，配送中心派出的送货车辆总数；
$f_v$：单次配送任务中，所有配送车辆行驶里程数之和；
$f_b$：单次配送任务中，所有配送线路上产生的货损量之和。

4.3.1 分析

决策变量：
$$X_{ijk}=\{\begin{array}{r}1,配送车k从市场(或配送中心)i到市场j\\ 0,其他\end{array}$$ $$Y_{ik}=\{\begin{array}{r}1,市场i的货物由配送车k配送\\ 0,其他\end{array}$$ 目标函数 $$minZ=min(f_0{c}_0+f_v{c}_v+f_{b}p)$$ $$\begin{gathered} minZ=min[k{c}_0+\sum _{i=0}^n\sum _{j=0}^n\sum _{k=1}^m{d}_{ij}{X}_{ijk}{c}_v \\ +(\sum _{i=0}^n\sum _{j=0}^n\sum _{k=1}^m{\beta }_1{t}_{ij}{u}_{ijk}{X}_{ijk}+\sum _{i=1}^n\sum _{k=1}^m{\beta }_{2}T{q}_i{Y}_{ik})p] \end{gathered}$$ 约束条件
每条路线的生鲜需求量之和不超过配送车最大装载量 $Q$：$$\sum _{i=1}^n{q}_i{Y}_{ik}\le Q(k=1,2,\cdots ,m)$$ 每辆车不可能同时去两个市场：$$\sum _{i=1}^n{X}_{oik}=\sum _{i=1}^n{X}_{iok}\le 1(k=1,2,\cdots ,m)$$ 每个市场只有一辆配送车送货：$$\sum _{k=1}^m{Y}_{ik}=1(i=1,2,\cdots ,n)$$ 假设各市场坐标及需求量如下表

序号	$x$	$y$	$q$
1	83	263	137
2	90	302	155
3	67	300	245
4	126	273	126
5	105	102	127
6	75	246	158
7	87	320	109
8	129	72	242
9	71	332	123
10	77	153	206
11	97	136	173
12	85	286	136
13	113	173	111
14	101	301	230
15	100	324	222

配送中心坐标为 $(93,245)$，配送车共7辆，单位里程费用 $1元/km$，每次配送固定成本为 $100元/辆$，平均行驶速度为 $50km/h$，最大载重量为 $500kg$，单位配送运输损耗比例为 $1\%$，单位装卸损耗比例为 $0.5\%$，生鲜装卸时间为 $15min$，单位重量货物的货损价格为 $10元/kg$；

编写程序，进行计算

4.4 TSP问题求解

这里用 4.3 的例子坐标进行展示，先编写所需要的 $matlab$ 函数：
① 让解在迭代中变化，相当于通过 $x_1$ 随机游走生成 $x_2$

% tsp_sa 随机扰动子函数
% 也就是原先的改变步长，这里通过交换两个地点的位置来进行计算
function circle = perturbation(circle)
    r = 1 + randperm(length(circle)-2);
    tmp = circle(r(1));
    circle(r(1)) = circle(r(2));
    circle(r(2)) = tmp;
end

② 计算路径总长度，函数输入值为距离矩阵 $d$ 和路径 $circle$

% tsp_sa 路线总长度子函数
function len = pathlength(d, circle)
    j = length(circle);
    len = 0;
    for i = 1:j-1
        len = len+d(circle(i), circle(i+1));
    end
end

③ 模拟退火主函数

% tsp_sa 模拟退火解决 TSP 问题
function [mincircle, minlen] = tsp_sa(d, cool, rho, iter)
    % d 是距离矩阵
    % rho 是衰减系数
    n = size(d, 1);                               % 地点数量
    temperature = 100 * n;                        % 初始温度
    circle = [1:n, 1];                            % 初始值
    mincircle = circle;                           % mincircle 是最小路径
    minlen = inf;                                 % minlen 是路径的长度
    s = 0;
    while temperature > cool                      % cool 是冷却温度
        for i= 1:iter                             % iter 是内循环次数
            len1 = pathlength(d, circle);         % 计算原路线总长度
            newcircle = perturbation(circle);     % 产生随机扰动
            len2 = pathlength(d, newcircle);      % 计算新路线总长度
    %---------------------------Metropolis------------------------------
            delta = len2 - len1;
            if exp(-delta/temperature) > rand
                circle = newcircle;
                len1 = len2;
            end
            if len1 < minlen
                mincircle = circle;
                minlen = len1;
            end
        end
        s = s + 1;
        temperature = temperature*rho;
    end
end

利用上述函数进行计算，得到结果如下：

clear;clc;
% ----------基础数据---------
market = [
    1	83	263	137;
    2	90	302	155;
    3	67	300	245;
    4	126	273	126;
    5	105	102	127;
    6	75	246	158;
    7	87	320	109;
    8	129	72	242;
    9	71	332	123;
    10	77	153	206;
    11	97	136	173;
    12	85	286	136;
    13	113	173	111;
    14	101	301	230;
    15	100	324	222];
center = [93 245];
A = market(:, 2:3);
B = [center; A];
C = num2str([0 1:size(A, 1)]');

h = pdist(B);
d = squareform(h);   % 距离矩阵
% --------------模拟退火初始控制变量---------------
cool = 0.1;
rho = 0.99;
iter = 1000;
% ----------------------求解---------------------
[mincircle, minlen] = tsp_sa(d, cool, rho, iter);

% ----------------------画图---------------------
scatter(B(:,1), B(:, 2), 'ro');
hold on;
plot(B(mincircle, 1), B(mincircle,2), 'g');
hold on;
text(B(:,1), B(:, 2), C);
title(['回路总长度', num2str(minlen)])

4.5 车辆路径优化

现有一大型零售企业的物流中心位于坐标 $(12,12)$(单位：$km$)，要给 $10$ 个超市送货，每个超市的需求量和坐标如下表所示。物流中心共 $6$ 辆车，每辆车载质量为 $40t$，要给这些超市送货，送货后回到起点。每辆车出行基本费用为 $200$ 元，单位运输成本为 $50元/km$，按直线距离计算路程。求配送费用最低的配送方案，用模拟退火算法求解。

超市编号	需求量	X	Y
1	10	3	13
2	15	12	7
3	21	2	26
4	9	24	12
5	25	31	22
6	18	5	14
7	23	17	11
8	15	16	13
9	8	12	8
10	13	8	2

（1）主函数

clc;clear;
%% 基础数据
center = [12, 12];                   % 物流中心
market = [1, 10, 3,  13;             % 超市的需求量和位置
          2, 15, 12,  7;
          3, 21, 2,  26;
          4, 9, 24,  12;
          5, 25, 31, 22;
          6, 18, 5,  14;
          7, 23 17,  11;
          8, 15, 16, 13;
          9, 8,  12,  8;
          10, 13, 8,  2;];
demands = market(:,2);               % 超市需求量  
cumsum = size(market, 1);            % 超市个数
market_xy = market(:, 3:4);          % 超市坐标

car = 6;                             % 最大车辆数
heavy = 40;                          % 每辆车的最大载重量
cost_0 = 200;                        % 每辆车出行基本费用
cost_v = 50;                         % 单位运输成本
punish = 3000;                       % 单辆车违反约束条件的惩罚

A = [center; market_xy];
h = pdist(A);
d = squareform(h);                   % 距离矩阵


% --------------模拟退火初始控制变量---------------
cool = 0.1;
rho = 0.99;
iter = 200;

% -----------------------------------------------------------
temperature = 100 * size(d, 1);                               % 初始温度

init_vc = init(cumsum, demands, heavy);                       % 路径初始化
[distance_all, n] = path_dist(init_vc, d);                    % 计算路径总距离和用车数量
mincircle = init_vc;
minlen = inf;                                                 % minlen 是路径的长度
% 初始化成功
vc = init_vc;

s = 0;
Y = [];
while temperature > cool                              % cool 是冷却温度
    for i= 1:iter                                     % iter 是内循环次数
        
        % 计算每个路径的距离以及路径的总距离
        [distance_all, n] = path_dist(init_vc, d);                    % 计算路径总距离和用车数量
        k = judge(vc, demands, heavy);                                % 判断是否有违约路径
        all_cost1 = cost_0 * n + cost_v * distance_all + k * punish;  % 计算成本
        path = get_path(vc, cumsum, car);                             % 将元胞拼接成路径
        
        % 获得新解 path2
        path2 = new_path(path);           
        new_vc = cut_path(path2, cumsum, car);                          % 将路径切割为元胞
        [distance_all2, n] = path_dist(new_vc, d);                      % 计算路径总距离和用车数量
        k2 = judge(new_vc, demands, heavy);
        all_cost2 = cost_0 * n + cost_v * distance_all2 + k2 * punish;  % 计算成本
        
%% ---------------------------Metropolis-------------------------------
        delta = all_cost2 - all_cost1;
        if exp(-delta/temperature) > rand
            vc = new_vc;
            all_cost1 = all_cost2;
        end
        if all_cost1 < minlen
            mincircle = vc;
            minlen = all_cost1;
        end
    end
    s = s + 1;
    Y(s) = minlen;
    temperature = temperature*rho;
end
minlen
mincircle
plot(Y, 'r')
ylabel('费用'); xlabel('迭代次数');

（2）路径初始化

% 初始化种群
function [init_vc] = init(cusnum, demands, heavy)
    j = ceil(rand*cusnum);
    k = 1;
    init_vc = cell(k, 1);
    if j == 1
        seq = 1:cusnum;
    elseif j == cusnum
        seq = [cusnum, 1:j-1];
    else
        seq1 = 1:j-1;
        seq2 = j:cusnum;
        seq = [seq2, seq1];
    end

    route = [];
    load = 0;
    i = 1;
    while i <= cusnum
        if load + demands(seq(i)) <= heavy           % 满足约束条件则加入路线，不满足则新开一条路线
            load = load + demands(seq(i));
            if isempty(route)
                route = [seq(i)];
            else
                route = [route, seq(i)];
            end
            if i == cusnum
                init_vc{k, 1} = route;
                break
            end
            i = i + 1;
        else
            init_vc{k, 1} = route;
            route = [];
            load = 0;
            k = k + 1;
        end
    end
end

（3）路径长度计算

% 计算路径的总长度
function [distance_all, n] = path_dist(vc, d)
    distance_all = zeros(1, length(vc));
    n = length(vc);
    for i = 1:n
        lr = vc{i, 1};
        dist = 0;
        for j = 1:length(lr)-1
            dist = dist + d(lr(j)+1, lr(j+1)+1);
        end
        dist = dist + d(1, lr(1)+1) + d(lr(end)+1, 1);
        distance_all(i) = dist;
    end
    distance_all = sum(distance_all);
end

（4）约束条件判断

% 判断是否违反约束条件
function k = judge(vc, demands, heavy)
    k = 0;
    for i = 1:length(vc)
        route = vc{i, 1};
        if sum(demands(route)) > heavy
            k = k + 1;
        end
    end
end

（5）获得新路径

%% 生成新路径
function path2 = new_path(path)
    n = length(path);
    path1 = path;
    c = randi(n, 1, 2);                  % 生成1~n中的两个随机整数用于交换
    p = 0.6;
    if rand < p                          % 交换法产生新路径，若随机概率落在目标概率内，则交换路径中的两个元素位置
        path2 = path1;                 
        path2(c(1)) = path1(c(2));
        path2(c(2)) = path1(c(1));
    else
        path2 = fliplr(path1);           % 左右翻转 
end

（6）切割一整条路径为多个子路径

% 路径切割
function vc = cut_path(path, cumsum, car)
    vc = cell(car, 1);
    n = find(path > cumsum);
    for i = 1:length(n)
       if i == 1
           vc{i,1} = path(1:n(i));
           vc{i,1}(end) = [];
       else 
           vc{i,1} = path(n(i-1):n(i));
           vc{i,1}(1) = [];
           vc{i,1}(end) = [];
       end
    end
    vc{end, 1} = path(n(end):end);
    vc{end, 1}(1) = [];
    vc(cellfun(@isempty, vc)) = [];
end

（7）将多个子路径拼接为一整个路径

% 路线拼接
function path = get_path(vc, cumsum, car)
    n = cumsum + car - 1;
    path = [];
    for i = 1:length(vc)
        rount = vc{i, 1};
        if i == length(vc)
            path = [path rount];
        else
            path = [path rount cumsum + i];
        end
    end
    if length(path) ~= n
        path(length(path)+1:n) = (length(path)+1):n; 
    end
end

结果得到：5辆车，-> 7 -> 8 -> ，-> 3 ->，-> 1 -> 6 ->，-> 5 -> 4 ->，-> 9 -> 2 -> 10 ->，花费为 $7614.5$ 元。

5.优缺点及可改进方向

5.1 优点

1、只需要单个初值就能进行计算，不像群体智能优化算法需要进行初始化种群；
2、可以找到全局最优点，类似梯度下降，但是梯度下降可能会被困在局部解，但是模拟退火有概率突跳性，虽然在温度较高时，接受差解的概率大，但当温度下降时，接受差解的概率就越来越小。

5.2 缺点

1、收敛速度慢，因为初始温度一般设定得很高，终止温度设定得低，这样才能体现模拟退火的特点，而且接受差解，并不是全程收敛，此时就需要时间跳出来；
2、受到温度、退火速度等模拟退火参数影响较大，比如衰减速度(退火速度)太大，可能会导致找不到全局最优解。

5.3 可改进方向

1、可以和其他智能算法结合在一起，比如模拟退火算法负责前半部分，遗传算法负责后半部分；又或者多种算法融合/混合在一起；
2、可以尝试多次找到比较好的模拟退火参数来作为最后使用的值。

6.参考文献

[1] 周佳莉. 模拟退火算法的应用[J]. 西部皮革,2019,41(20):82.
[2] 包子阳,余继周,杨杉. 智能优化算法及其MATLAB实例（第3版）
[3] 康雯轩. 基于模拟退火算法的共享单车城市配送路径规划[J]. 科技与创新,2022(13):104-106,109. DOI:10.15913/j.cnki.kjycx.2022.13.032.
[4] 胡良剑,孙晓君. MATLAB数学实验（第3版）