模拟退火算法学习笔记
模拟退火算法(SA):是一种模拟物理退火过程而设计的优化算法。先在一个高温状态下,然后逐渐退火,在每个温度下慢慢冷却,最终达到物理基态(相当于算法找到最优解)
作用;求解TSP问题、求最值、全局优化、生产调度、控制工程、机器学习、信号处理等问题
爬山算法(贪心算法):对问题求解时,总是做出当前看起来是最好的选择,,不考虑整体最优,得到的一般是局部最优解。
模拟退火算法:也是贪心算法,但是在其过程中引入了随机因素,以一定概率接受一个比当前解要差的解。
模拟退火算法源于对固体退火过程的模拟,采用Metropolis准则,并用一组称为冷却进度表的参数控制算法的进程,使得算法在多项式时间里可以给出一个近似最优解。
模拟退火算法在某一初温下,伴随温度的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解,记载局部最优解的空间内能概率性地跳出并最终趋于全局最优。
物理退火过程:
- 在初始状态,固体内部的粒子存在非均匀状态
- 加热到一定温度,增强其粒子的热运动,使粒子偏离其平衡位置,消除原来的非均匀状态
- 让温度慢慢降低,达到热平衡状态,粒子逐渐均匀有序,最终达到内能最小的状态
模拟退火过程:
- 设定当前解(即为当前的最优解):令
,即开始退火的初始温度,随机生成一个初始解
,并计算相应的目标函数值
。
- 产生新解与当晚解差值:根据当前解
进行扰动,产生一个新解
,计算相应的目标函数值
,得到
。
- 判断新解是否被接受:若
,则新解
被接受,若
,则新解
按概率
接受,
为当前温度。(一般模拟退火算法求解最小值)
- 当新解被确定接受时:新解
被当作当前解
- 循环以上四个步骤:在温度
下,重复
次的扰动和接受过程,接着执行下一步骤。
- 最后找到全局最优解:判断
是否已经达到终止温度
,是,则终止算法;否,则跳转到循环步骤继续执行。
冷却进度表:退火过程由一组初始参数,即冷却进度表 控制。他的目的是尽量使系统达到平衡,以使算法在有限时间内逼近最优解,包括:
- 控制温度参数的初值
- 控制温度
的衰减函数(温度的更新)
- 马尔科夫链的长度
(迭代次数)
- 控制参数
的终值(停止准则)
实验表明:初温越大,获得全局最优解的几率越大,但花费的时长越长
Metropolis准则:在某个温度下固体分子从一个状态转移到另一个状态时,如果新状态的内能小,则无条件接受;如果新状态的内能大,则以一定的概率接受它。
马尔科夫链的长度;任意温度
的迭代次数 算法在马尔科夫链长度内持续进行“产生新解一判断一接受/舍弃”的迭代过程,对应着固体在某-恒定温度 下趋于热平衡的过程。若在一定的温度下做无限次迭代,相应的马尔科夫链可以达到平稳分布概率。 马尔科夫链的选取还与温度控制参数
的下降密切相关,缓慢下降可以避免过长的马尔科夫链。在控制参数的衰减函数已经选定的前提下,让每个取值都能够达到准平衡状态。 根据这一原则-般取
, 其中
为问题的规模。
控制参数的终值
(停止准则): 合理的停止准则既能保证算法收敛于某一近似解,又能使最终解具有一定的全局性。通常可以根据迭代的次数或者终止的温度或者迭代过程在若干个相继的链中的解没有任何变化等条件来判断迭代的终止。
最终温度通常是0,但会耗掉许多模拟时间。温度趋近于0,其周围状态几乎是一样的。所以找寻一个低到可接受的温度即可。
控制温度的衰减函数(温度的更新): 不同退火方法的温度下降速度是不一样的,其中指数降温是最常用的一种退火策略,其温度变化很有规律,直接与参数相关,使我们研究的主要对象。
常用的衰减函数:
其中是一个接近1的常数。一般取
。该衰减函数对控制参数的衰减量是随算法进程递减的。
注意:小的衰减量可能导致算法进程迭代次数的增加,从而使算法进程接受更多的变换,从而访问更多的临域,搜索更大范围的解空间,返回更高质量的最终解。
例题一:TSP(旅行商问题)
一位旅行者从一出发点出发,要求经过31个城市(目标点的坐标已给出),并且每个点只能经过一次,最终经过所有点后返回起点。要求:为旅行者制定一条最短路径。
%%%%%%%%%%%%%%%%%%%%%%模拟退火算法解决TSP问题%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
clear all; %清除所有变量
close all; %清图
clc; %清屏
C=[1304 2312; % 城市坐标
3639 1315;
4177 2244;
3712 1399;
3488 1535;
3326 1556;
3238 1229;
4196 1044;
4312 790;
4386 570;
3007 1970;
2562 1756;
2788 1491;
2381 1676;
1332 695;
3715 1678;
3918 2179;
4061 2370;
3780 2212;
3676 2578;
4029 2838;
4263 2931;
3429 1980;
3507 2376;
3394 2643;
3439 3201;
2935 3240;
3140 3550;
2545 2357;
2778 2826;
2370 2975]; %31个城市坐标
n=size(C,1); %TSP问题的规模,即城市数目
T=1000; %初始温度
L=100; %马可夫链长度
K=0.99; %衰减参数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%城市坐标结构体%%%%%%%%%%%%%%%%%%%%%%%%%%
city=struct([]); %结构体变量,类似python中的字典
for i=1:n %city(i)的值为第i座城市的坐标
city(i).x=C(i,1);
city(i).y=C(i,2);
end
l=1; %统计迭代次数
len(l)=func3(city,n); %每次迭代后的路线长度
% figure(1);
while T>0.001 %终止温度
%%%%%%%%%%%%%%%%多次迭代扰动,温度降低之前多次实验%%%%%%%%%%%%%%%
for i=1:L
%%%%%%%%%%%%%%%%%%%计算原路线总距离%%%%%%%%%%%%%%%%%%%%%%%%%
len1=func3(city,n);
%%%%%%%%%%%%%%%%%%%%%%%%%产生随机扰动%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%随机置换两个不同的城市的坐标%%%%%%%%%%%%%%%%%
p1=floor(1+n*rand()); %朝负无穷方向取整,如floor(-1.3)=-2,ceil相反
p2=floor(1+n*rand()); %
while p1==p2
p1=floor(1+n*rand());
p2=floor(1+n*rand());
end
tmp_city=city;
tmp=tmp_city(p1);
tmp_city(p1)=tmp_city(p2);
tmp_city(p2)=tmp;
%%%%%%%%%%%%%%%%%%%%%%%%计算新路线总距离%%%%%%%%%%%%%%%%%%%%
len2=func3(tmp_city,n);
%%%%%%%%%%%%%%%%%%新老距离的差值,相当于能量%%%%%%%%%%%%%%%%%
delta_e=len2-len1;
%%%%%%%%%%%%新路线好于旧路线,用新路线代替旧路线%%%%%%%%%%%%%%
if delta_e<0
city=tmp_city;
else
%%%%%%%%%%%%%%%%%%以概率选择是否接受新解%%%%%%%%%%%%%%%%%
if exp(-delta_e/T)>rand()
city=tmp_city;
end
end
end
l=l+1;
%%%%%%%%%%%%%%%%%%%%%%%%%计算新路线距离%%%%%%%%%%%%%%%%%%%%%%%%%%
len(l)=func3(city,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%温度不断下降%%%%%%%%%%%%%%%%%%%%%%%%%%
T=T*K;
% for i=1:n-1
% plot([city(i).x,city(i+1).x],[city(i).y,city(i+1).y],'bo-');
% hold on;
% end
% plot([city(n).x,city(1).x],[city(n).y,city(1).y],'ro-');
% title(['优化最短距离:',num2str(len(l))]);
% hold off;
% pause(0.005);
end
figure(1);
for i=1:n-1
plot([city(i).x,city(i+1).x],[city(i).y,city(i+1).y],'bo-');
hold on;
end
plot([city(n).x,city(1).x],[city(n).y,city(1).y],'ro-');
title(['优化最短距离:',num2str(len(l))])
hold off;
figure(2);
plot(len)
% plot(C(:,1),C(:,2),'bo-')以C的第1列为横坐标,第二列为纵坐标
xlabel('迭代次数')
ylabel('目标函数值')
title('路径长变化曲线')
toc
function len=func3(city,n)
len=0;
for i=1:n-1
len=len+sqrt((city(i).x-city(i+1).x)^2+(city(i).y-city(i+1).y)^2);
end
%len=len+sqrt((city(1).x-start(1))^2+(city(1).y-start(2))^2);
len=len+sqrt((city(n).x-city(1).x)^2+(city(n).y-city(1).y)^2);
end
最终结果(近似结果,每次运行可能结果不同)
例题二:
求解一元函数的最小值,自变量x的范围[1,2]。
使用MATLAB工具箱(APP):
function costVal=cost(x)
costVal=sin(5*pi*x)/x;
end
最后结果:
例题三:指派问题
每行代表A、B、C、D、E、F每个人单独做甲、乙、丙、丁、戊、己六个工作需要的花费。每列代表每个工作分别由每个人完成需要的花费。 即:矩阵中位于(i,j)的元素是第i个人做第j个工作的花费。试问怎么安排最合理?
% MarkoveLength 马尔科夫链长度
% DecayScale 温度衰减参数
clc
clear all
MarkoveLength=1000; %迭代次数
DecayScale=0.9; % 衰减参数α
Temperature=1000; % initial temperature
t = 1; % final temperature
% 指派矩阵,每个人做每件事所需要的时间
assingnMatrix=[17 14 2 10 1 20;
5 8 4 18 12 9;
1 12 2 10 1 20;
20 6 15 9 7 11;
17 15 13 10 10 20
12 7 20 12 9 10];
n=length(assingnMatrix);
%% 初始解
x=[1,2,3,4,5,6];
totalCost = 0;
for i=1:n
totalCost=totalCost+assingnMatrix(i,x(i));%计算总的时间
end
BestCost=totalCost;
BestAssign=x;
while (Temperature > t || abs(deltaCost)<2)
for i = 1:MarkoveLength
r = randperm(6,2); %产生两个随机数,用来交换x中的任务分配顺序
r1=r(1);
r2=r(2);
temp=x(r1);
x(r1)=x(r2);
x(r2)=temp;
% 计算
totalCost=0;
for k=1:n
totalCost=totalCost+assingnMatrix(k,x(k));
end
% 判断时间是否优化
deltaCost=totalCost-BestCost;%f(x')-f(x)
if deltaCost <= 0
BestCost = totalCost; % 若时间减少则无条件接受
BestAssign=x;
else
if (rand()<exp(-deltaCost/Temperature))
BestCost = totalCost; %否则在一定概率下接受
BestAssign=x;
end
end
end
Temperature=Temperature*DecayScale;
end
disp(BestCost);
disp(assingnMatrix);
disp(BestAssign);