matlab智能算法之模拟退火算法

智能算法之模拟退火算法

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 组合优化与物理退化

3.原理

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

3.1 算法核心迭代

设 fk,fk+1f_k,f_{k+1}fk​,fk+1​ 分别为目标函数在第 kkk 次和第 k+1k+1k+1 次迭代值，即 fk=f(Xk),fk+1=f(Xk+1)f_k=f(X_k),f_{k+1}=f(X_{k+1})fk​=f(Xk​),fk+1​=f(Xk+1​)。

若 fk>fk+1f_k>f_{k+1}fk​>fk+1​，则接受 Xk+1X_{k+1}Xk+1​ 作为下一次迭代的初值继续进行运算，直至满足终止条件（收敛或者迭代次数耗尽）；

若 fk<fk+1f_k<f_{k+1}fk​<fk+1​，则接受 Xk+1X_{k+1}Xk+1​ 的概率为 ppp，拒绝的概率为 1−p1-p1−p，其中接受概率公式如下：p=exp(−fk+1−fkT)p=exp(-\frac{f_{k+1}-f_k}{T})p=exp(−Tfk+1​−fk​​) 其中 TTT 为控制参数，在模拟退火算法中，TTT 必须缓慢减少，避免变化过快，导致优化陷入极值点。

3.2 具体流程

以求解极小值问题为例：

（1）（1）（1）设定初始温度 TTT，迭代次数 LLL，误差 ϵ\epsilonϵ，初始解 X0X_0X0​，得到目标函数f(X0)f(X_0)f(X0​)；

（2）（2）（2）循环进行迭代，k=1,⋯,Lk=1,\cdots,Lk=1,⋯,L；

（3）（3）（3）根据当前温度得到新解 X′X'X′(第一个可以随机生成)， 代入 f(X)f(X)f(X) 得到 f(X′)f(X')f(X′)；

（4）（4）（4）若 f(X)>f(X′)f(X)>f(X')f(X)>f(X′)，则接受 X′X'X′ 作为下一次迭代的初值，若 f(X)<f(X′)f(X)<f(X')f(X)<f(X′)，则以概率 exp(−f(X′)−f(X)T)exp(-\frac{f(X')-f(X)}{T})exp(−Tf(X′)−f(X)​) 接受 X′X'X′ 作为下一次迭代的初值；

（5）（5）（5）若满足终止条件，则输出当前最优解 X′X'X′，终止条件；

（6）（6）（6）TTT 减小，然后回到流程（2）（2）（2），在当前温度下继续循环迭代，执行步骤（3，4，5）（3，4，5）（3，4，5）；

（7）（7）（7）注意，这里的终止条件可以有多个，比如：迭代次数耗尽、满足最优解、温度降低到一定值、新目标函数与原先目标函数的误差小于 ϵ\epsilonϵ。

4.案例

4.1 求解n元函数的极小值

例1：求解函数 f(X)=X2f(X)=X^2f(X)=X2 的最小值，其中 X2=∑i=110x2X^2=\sum_{i=1}^{10}x^2X2=∑i=110​x2，−10≤x≤10-10\leq x\leq 10−10≤x≤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，极大值问题为-1n = 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:Lx2 = x1 + s*(rand(n, 1)*(up-sup)-up);%-----下面这个循环是写题目的条件的，这里没啥好些，就重复了边界条件for j = 1:nwhile (x2(j) > up) || (x2(j) < sup)x2(j) = x1(j) + s*(rand*(up-sup)-up);endend%--------------查看解是否是全局最优解----------------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;elsep = exp(-(func1(x2)-func1(x1))/T0);if p > rand(1)x1 = x2;P = P + 1;endendY(P) = func1(NowBest);enddelta = 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]X=[-0.0114,-0.0002,0.0020,0.0101,0.0074,-0.0003,-0.0084,0.0112,-0.0252,0.0052]X=[−0.0114,−0.0002,0.0020,0.0101,0.0074,−0.0003,−0.0084,0.0112,−0.0252,0.0052]

f(X)=0.0011f(X)=0.0011f(X)=0.0011

绘制适应度进化曲线：

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

4.2 求解二元函数的极小值

例2：求解函数 f(x,y)=5cos(xy)+xy+y3f(x,y) = 5cos(xy)+xy+y^3f(x,y)=5cos(xy)+xy+y3 的最小值，其中 −5≤x≤5,−5≤y≤5-5\leq x\leq 5,-5\leq y\leq 5−5≤x≤5,−5≤y≤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)f(x,y)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:Nfor j = 1:Nz(i, j)=5*cos(x1(i)*y1(j))+x1(i)*y1(j)+y1(j)^3;endendmesh(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,y0func = @(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.0788f(x,y)=f(4.4242,-5.0000)=-152.0788f(x,y)=f(4.4242,−5.0000)=−152.0788

4.3 货物配送规划

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

nnn ：销售市场数量；

qiq_iqi​：每个市场需求量 (i=1,2,⋯,n)(i=1,2,\cdots,n)(i=1,2,⋯,n)；

mmm：配送车数量（每辆车都型号都一致）；

QQQ：配送车最大装载量；

dijd_{ij}dij​：市场 iii 到市场 jjj 的距离；

tijt_{ij}tij​：从市场 iii 到市场 jjj 所需时间；

doid_{oi}doi​：配送中心 ooo 到市场 iii 的距离；

c0c_0c0​ ：每次配送固定成本；

cvc_vcv​：单位里程行驶费用；

vvv：配送车匀速行驶速度；

TTT：生鲜装卸时间；

β1\beta_1β1​：单位时间配送运输损耗比例；

β2\beta_2β2​：单位时间配送装卸损耗比例；

uijku_{ijk}uijk​：配送车辆 kkk 在 iii 到 jjj 路径上行驶时的载重量；

ppp：单位重量货物的货损价格；

f0f_0f0​：单次配送任务中，配送中心派出的送货车辆总数；

fvf_vfv​：单次配送任务中，所有配送车辆行驶里程数之和；

fbf_bfb​：单次配送任务中，所有配送线路上产生的货损量之和。

4.3.1 分析

决策变量：

Xijk={1,配送车k从市场(或配送中心)i到市场j0,其他X_{ijk}=\left\{ \begin{aligned} 1, \quad 配送车 k 从市场(或配送中心) i 到市场 j \\ 0, \quad \quad \quad \quad \quad \quad \quad \quad其他 \quad \quad \quad\quad \quad \quad \quad \end{aligned} \right. Xijk​={1,配送车k从市场(或配送中心)i到市场j0,其他​ Yik={1,市场i的货物由配送车k配送0,其他Y_{ik}=\left\{ \begin{aligned} 1, \quad 市场 i 的货物由配送车k配送 \\ 0, \quad \quad \quad \quad \quad其他\quad \quad \quad \quad\quad \quad \end{aligned} \right. Yik​={1,市场i的货物由配送车k配送0,其他​目标函数minZ=min(f0c0+fvcv+fbp)min Z = min(f_0c_0+f_vc_v+f_bp)minZ=min(f0​c0​+fv​cv​+fb​p) minZ=min[kc0+∑i=0n∑j=0n∑k=1mdijXijkcv+(∑i=0n∑j=0n∑k=1mβ1tijuijkXijk+∑i=1n∑k=1mβ2TqiYik)p]minZ=min[kc_0+\sum_{i=0}^n\sum_{j=0}^n\sum_{k=1}^md_{ij}X_{ijk}c_v \\ + (\sum_{i=0}^n\sum_{j=0}^n\sum_{k=1}^m\beta_1t_{ij}u_{ijk}X_{ijk}+\sum_{i=1}^n\sum_{k=1}^m\beta_2Tq_{i}Y_{ik})p]minZ=min[kc0​+i=0∑n​j=0∑n​k=1∑m​dij​Xijk​cv​+(i=0∑n​j=0∑n​k=1∑m​β1​tij​uijk​Xijk​+i=1∑n​k=1∑m​β2​Tqi​Yik​)p]约束条件

每条路线的生鲜需求量之和不超过配送车最大装载量 QQQ：∑i=1nqiYik≤Q(k=1,2,⋯,m)\sum_{i=1}^nq_{i}Y_{ik} \leq Q \quad(k=1,2,\cdots,m)i=1∑n​qi​Yik​≤Q(k=1,2,⋯,m) 每辆车不可能同时去两个市场：∑i=1nXoik=∑i=1nXiok≤1(k=1,2,⋯,m)\sum_{i=1}^nX_{oik}=\sum_{i=1}^nX_{iok}\leq1 \quad(k=1,2,\cdots,m)i=1∑n​Xoik​=i=1∑n​Xiok​≤1(k=1,2,⋯,m) 每个市场只有一辆配送车送货：∑k=1mYik=1(i=1,2,⋯,n)\sum_{k=1}^mY_{ik}=1\quad (i=1,2,\cdots,n)k=1∑m​Yik​=1(i=1,2,⋯,n) 假设各市场坐标及需求量如下表

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

编写程序，进行计算

4.4 TSP问题求解

这里用4.3的例子坐标进行展示，先编写所需要的 matlabmatlabmatlab 函数：

① 让解在迭代中变化，相当于通过 x1x_1x1​ 随机游走生成 x2x_2x2​

% 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

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

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

③ 模拟退火主函数

% 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) > randcircle = newcircle;len1 = len2;endif len1 < minlenmincircle = circle;minlen = len1;endends = s + 1;temperature = temperature*rho;endend

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

clear;clc;% ----------基础数据---------market = [183263137;290302155;367300245;4126273126;5105102127;675246158;78739;812972242;971332123;1077153206;1197136173;1285286136;13113173111;14101301230;15100324222];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)(12,12)(12,12)(单位：kmkmkm)，要给 101010 个超市送货，每个超市的需求量和坐标如下表所示。物流中心共 666 辆车，每辆车载质量为 40t40t40t，要给这些超市送货，送货后回到起点。每辆车出行基本费用为 200200200 元，单位运输成本为 50元/km50元/km50元/km，按直线距离计算路程。求配送费用最低的配送方案，用模拟退火算法求解。

（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); % 将元胞拼接成路径% 获得新解 path2path2 = 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) > randvc = new_vc;all_cost1 = all_cost2;endif all_cost1 < minlenmincircle = vc;minlen = all_cost1;endends = s + 1;Y(s) = minlen;temperature = temperature*rho;endminlenmincircleplot(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 == 1seq = 1:cusnum;elseif j == cusnumseq = [cusnum, 1:j-1];elseseq1 = 1:j-1;seq2 = j:cusnum;seq = [seq2, seq1];endroute = [];load = 0;i = 1;while i <= cusnumif load + demands(seq(i)) <= heavy % 满足约束条件则加入路线，不满足则新开一条路线load = load + demands(seq(i));if isempty(route)route = [seq(i)];elseroute = [route, seq(i)];endif i == cusnuminit_vc{k, 1} = route;breakendi = i + 1;elseinit_vc{k, 1} = route;route = [];load = 0;k = k + 1;endendend

（3）路径长度计算

% 计算路径的总长度function [distance_all, n] = path_dist(vc, d)distance_all = zeros(1, length(vc));n = length(vc);for i = 1:nlr = vc{i, 1};dist = 0;for j = 1:length(lr)-1dist = dist + d(lr(j)+1, lr(j+1)+1);enddist = dist + d(1, lr(1)+1) + d(lr(end)+1, 1);distance_all(i) = dist;enddistance_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)) > heavyk = k + 1;endendend

（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));elsepath2 = 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 == 1vc{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) = [];endendvc{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];elsepath = [path rount cumsum + i];endendif length(path) ~= npath(length(path)+1:n) = (length(path)+1):n; endend

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

5.优缺点及可改进方向

5.1 优点

1、只需要单个初值就能进行计算，不像群体智能优化算法需要进行初始化种群；

2、可以找到全局最优点，类似梯度下降，但是梯度下降可能会被困在局部解，但是模拟退火有概率突跳性，虽然在温度较高时，接受差解的概率大，但当温度下降时，接受差解的概率就越来越小。

5.2 缺点

1、收敛速度慢，因为初始温度一般设定得很高，终止温度设定得低，这样才能体现模拟退火的特点，而且接受差解，并不是全程收敛，此时就需要时间跳出来；

2、受到温度、退火速度等模拟退火参数影响较大，比如衰减速度(退火速度)太大，可能会导致找不到全局最优解。

5.3 可改进方向

1、可以和其他智能算法结合在一起，比如模拟退火算法负责前半部分，遗传算法负责后半部分；又或者多种算法融合/混合在一起；

2、可以尝试多次找到比较好的模拟退火参数来作为最后使用的值。

6.参考文献

[1] 周佳莉. 模拟退火算法的应用[J]. 西部皮革,,41(20):82.

[2] 包子阳,余继周,杨杉. 智能优化算法及其MATLAB实例（第3版）

[3] 康雯轩. 基于模拟退火算法的共享单车城市配送路径规划[J]. 科技与创新,(13):104-106,109. DOI:10.15913/ki.kjycx..13.032.

[4] 胡良剑,孙晓君. MATLAB数学实验（第3版）