参考链接:模拟退火算法(SA算法)
【 1. 自然现象中的退火 】
-
退火
将固体加热到足够高的温度,使分子呈随机排列状态,然后逐步降温使之冷却,最后分子以低能状态排列,固体达到某种稳定状态。 -
物理退火过程
加温过程:增强粒子的热运动,消除系统原先可能存在的非均匀态;
等温过程:对于与环境换热而温度不变的封闭系统,系统状态的自发变化总是朝自由能减少的方向进行,当自由能达到最小时,系统达到平衡态;
冷却过程:使粒子热运动减弱并渐趋有序,系统能量逐渐下降,从而得到低能的晶体结构。 -
热力学中的退火现象
指物体逐渐降温时发生的物理现象:
温度越低,物体的能量状态越低,到达足够的低点时,液体开始冷凝与结晶,在结晶状态时,系统的能量状态最低。缓慢降温时,可达到最低能量状态;但如果快速降温,会导致不是最低能态的非晶形。 -
大自然现象
缓缓降温,使得物体分子在每一温度时,能够有足够时间找到安顿位置,则逐渐地,到最后可得到最低能态,系统最稳定。
【 2. 模拟退火算法思想 】
模仿自然界退火现象而得,利用了物理中固体物质的退火过程与一般优化问题的相似性。模拟退火算法从某一较高初温出发,伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解,即在局部最优解能概率性地跳出并最终趋于全局最优。
设热力学系统 S 中有 n 个状态(有限个、离散的 ),状态 i 的能量为
E
i
E_i
Ei,在温度
T
k
T_k
Tk 下,经一段时间达到热平衡,这时处于状态 i 的概率为:
P
i
(
T
k
)
=
C
k
×
e
−
E
i
T
k
=
C
k
×
e
1
E
i
T
k
P_i(T_k)=C_k\times e^{\frac{-E_i}{T_k}}=C_k\times e^{\frac{1}{\frac{E_i}{T_k}}}
Pi(Tk)=Ck×eTk−Ei=Ck×eTkEi1
则有如下关系:
E
i
↓
→
P
i
↑
E_i↓ → P_i↑
Ei↓→Pi↑
T
k
↓
→
P
i
↓
T_k↓ → P_i↓
Tk↓→Pi↓
Bolzman方程(描述非热力学平衡状态的热力学系统统计行为):
P
i
(
T
k
)
=
e
−
E
i
T
k
∑
i
=
1
n
e
−
E
i
T
k
P_i(T_k)=\frac{e^{\frac{-E_i}{T_k}}}{\sum_{i=1}^n{e^{\frac{-E_i}{T_k}}}}
Pi(Tk)=∑i=1neTk−EieTk−Ei
则有如下关系:
E
i
↓
→
P
i
(
T
k
)
↑
E_i↓ → P_i(T_k)↑
Ei↓→Pi(Tk)↑
T
k
↓
→
P
i
∗
(
T
k
)
↑
并
趋
于
1
,
i
∗
代
表
E
i
最
小
时
T_k↓→ P_{i^*}(T_k)↑并趋于1,i^*代表E_i最小时
Tk↓→Pi∗(Tk)↑并趋于1,i∗代表Ei最小时
平
均
能
量
方
程
:
E
‾
(
T
k
)
=
∑
i
=
1
n
E
i
⋅
P
i
(
T
k
)
平均能量方程:\overline{\text{E}}(T_k)=\sum_{i=1}^nE_i·P_i(T_k)
平均能量方程:E(Tk)=∑i=1nEi⋅Pi(Tk),
T
k
↓
→
E
‾
趋
于
E
i
∗
T_k↓→ \overline{\text{E}}趋于E_{i^*}
Tk↓→E趋于Ei∗
①
当
T
k
很
大
时
,
E
i
T
k
趋
于
0
,
P
i
(
T
k
)
=
C
k
=
定
值
=
1
n
,
此
时
各
状
态
概
率
相
等
。
① 当T_k很大时,\frac{E_i}{T_k}趋于0,P_i(T_k)=C_k=定值=\frac{1}{n},此时各状态概率相等。
①当Tk很大时,TkEi趋于0,Pi(Tk)=Ck=定值=n1,此时各状态概率相等。
②
当
T
k
很
小
趋
于
0
时
,
E
i
T
k
趋
于
∞
,
以
概
率
1
并
趋
于
最
小
能
量
状
态
。
② 当T_k很小趋于0时,\frac{E_i}{T_k}趋于∞,以概率1并趋于最小能量状态。
②当Tk很小趋于0时,TkEi趋于∞,以概率1并趋于最小能量状态。
- 模拟退火的基本思想
(1) 初始化:初始温度T(充分大),初始解状态S(是算法迭代的起点),每个T值的迭代次数L
(2) 对k=1, …, L做第(3)至第6步:
(3) 产生新解S′
(4) 计算增量ΔT=C(S′)-C(S),其中C(S)为评价函数
(5) 若ΔT<0则接受S′作为新的当前解,否则以概率exp(-ΔT/T)接受S′作为新的当前解.
(6) 如果满足终止条件则输出当前解作为最优解,结束程序。
终止条件通常取为连续若干个新解都没有被接受时终止算法。
(7) T逐渐减少,且T->0,然后转第2步。
邻域是一个特殊的区间,以点a为中心点任何开区间称为点a的邻域,记作U(a)。
【 3. 模拟退火算法对TSP问题的求解 】
旅行商问题,即TSP问题(Travelling Salesman Problem)又译为旅行推销员问题、货郎担问题,是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市,他必须选择所要走的路径,路径的限制是每个城市只能拜访一次,而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。
迄今为止,这类问题中没有一个找到有效算法。倾向于接受NP完全问题(NP-Complet或NPC)和NP难题(NP-Hard或NPH)不存在有效算法这一猜想,认为这类问题的大型实例不能用精确算法求解,必须寻求这类问题的有效的近似算法。
TSP“旅行商问题”的应用领域包括:如何规划最合理高效的道路交通,以减少拥堵;如何更好地规划物流,以减少运营成本;在互联网环境中如何更好地设置节点,以更好地让信息流动等。
【 4. 模拟退火算法的 Matlab描述 】
swap.m
function [ newpath , position ] = swap( oldpath , number )
% 对 oldpath 进 行 互 换 操 作
% number 为 产 生 的 新 路 径 的 个 数
% position 为 对 应 newpath 互 换 的 位 置
m = length( oldpath ) ; % 城 市 的 个 数
newpath = zeros( number , m ) ;
position = sort( randi( m , number , 2 ) , 2 ); % 随 机 产 生 交 换 的 位 置
for i = 1 : number
newpath( i , : ) = oldpath ;
% 交 换 路 径 中 选 中 的 城 市
newpath( i , position( i , 1 ) ) = oldpath( position( i , 2 ) ) ;
newpath( i , position( i , 2 ) ) = oldpath( position( i , 1 ) ) ;
end
pathfare.m
function [ objval ] = pathfare( fare , path )
% 计 算 路 径 path 的 代 价 objval
% path 为 1 到 n 的 排 列 ,代 表 城 市 的 访 问 顺 序 ;
% fare 为 代 价 矩 阵 , 且 为 方 阵 。
[ m , n ] = size( path ) ;
objval = zeros( 1 , m ) ;
for i = 1 : m
for j = 2 : n
objval( i ) = objval( i ) + fare( path( i , j - 1 ) , path( i , j ) ) ;
end
objval( i ) = objval( i ) + fare( path( i , n ) , path( i , 1 ) ) ;
end
distance.m
function [ fare ] = distance( coord )
% 根 据 各 城 市 的 距 离 坐 标 求 相 互 之 间 的 距 离
% fare 为 各 城 市 的 距 离 , coord 为 各 城 市 的 坐 标
[ v , m ] = size( coord ) ; % m 为 城 市 的 个 数
fare = zeros( m ) ;
for i = 1 : m % 外 层 为 行
for j = i : m % 内 层 为 列
fare( i , j ) = ( sum( ( coord( : , i ) - coord( : , j ) ) .^ 2 ) ) ^ 0.5 ;
fare( j , i ) = fare( i , j ) ; % 距 离 矩 阵 对 称
end
end
myplot.m
function [ ] = myplot( path , coord , pathfar )
% 做 出 路 径 的 图 形
% path 为 要 做 图 的 路 径 ,coord 为 各 个 城 市 的 坐 标
% pathfar 为 路 径 path 对 应 的 费 用
len = length( path ) ;
clf ;
hold on ;
title( [ '近似最短路径如下,路程为' , num2str( pathfar ) ] ) ;
plot( coord( 1 , : ) , coord( 2 , : ) , 'ok');
pause( 0.4 ) ;
for ii = 2 : len
plot( coord( 1 , path( [ ii - 1 , ii ] ) ) , coord( 2 , path( [ ii - 1 , ii ] ) ) , '-b');
x = sum( coord( 1 , path( [ ii - 1 , ii ] ) ) ) / 2 ;
y = sum( coord( 2 , path( [ ii - 1 , ii ] ) ) ) / 2 ;
text( x , y , [ '(' , num2str( ii - 1 ) , ')' ] ) ;
pause( 0.4 ) ;
end
plot( coord( 1 , path( [ 1 , len ] ) ) , coord( 2 , path( [ 1 , len ] ) ) , '-b' ) ;
x = sum( coord( 1 , path( [ 1 , len ] ) ) ) / 2 ;
y = sum( coord( 2 , path( [ 1 , len ] ) ) ) / 2 ;
text( x , y , [ '(' , num2str( len ) , ')' ] ) ;
pause( 0.4 ) ;
hold off ;
clear;
% 程 序 参 数 设 定
Coord = ... % 城 市 的 坐 标 Coordinates
[ 0.6683 0.6195 0.4 0.2439 0.1707 0.2293 0.5171 0.8732 0.6878 0.8488 ; ...
0.2536 0.2634 0.4439 0.1463 0.2293 0.761 0.9414 0.6536 0.5219 0.3609 ] ;
t0 = 1 ; % 初 温 t0
iLk = 20 ; % 内 循 环 最 大 迭 代 次 数 iLk
oLk = 50 ; % 外 循 环 最 大 迭 代 次 数 oLk
lam = 0.95 ; % λ lambda
istd = 0.001 ; % 若 内 循 环 函 数 值 方 差 小 于 istd 则 停 止
ostd = 0.001 ; % 若 外 循 环 函 数 值 方 差 小 于 ostd 则 停 止
ilen = 5 ; % 内 循 环 保 存 的 目 标 函 数 值 个 数
olen = 5 ; % 外 循 环 保 存 的 目 标 函 数 值 个 数
% 程 序 主 体
m = length( Coord ) ; % 城 市 的 个 数 m
fare = distance( Coord ) ; % 路 径 费 用 fare
path = 1 : m ; % 初 始 路 径 path
pathfar = pathfare( fare , path ) ; % 路 径 费 用 path fare
ores = zeros( 1 , olen ) ; % 外 循 环 保 存 的 目 标 函 数 值
e0 = pathfar ; % 能 量 初 值 e0
t = t0 ; % 温 度 t
for out = 1 : oLk % 外 循 环 模 拟 退 火 过 程
ires = zeros( 1 , ilen ) ; % 内 循 环 保 存 的 目 标 函 数 值
for in = 1 : iLk % 内 循 环 模 拟 热 平 衡 过 程
[ newpath , v ] = swap( path , 1 ) ; % 产 生 新 状 态
e1 = pathfare( fare , newpath ) ; % 新 状 态 能 量
% Metropolis 抽 样 稳 定 准 则
r = min( 1 , exp( - ( e1 - e0 ) / t ) ) ;
if rand < r
path = newpath ; % 更 新 最 佳 状 态
e0 = e1 ;
end
ires = [ ires( 2 : end ) e0 ] ; % 保 存 新 状 态 能 量
% 内 循 环 终 止 准 则 :连 续 ilen 个 状 态 能 量 波 动 小 于 istd
if std( ires , 1 ) < istd
break ;
end
end
ores = [ ores( 2 : end ) e0 ] ; % 保 存 新 状 态 能 量
% 外 循 环 终 止 准 则 :连 续 olen 个 状 态 能 量 波 动 小 于 ostd
if std( ores , 1 ) < ostd
break ;
end
t = lam * t ;
end
pathfar = e0 ;
% 输 入 结 果
fprintf( '近似最优路径为:\n ' )
%disp( char( [ path , path(1) ] + 64 ) ) ;
disp(path)
fprintf( '近似最优路径路程\tpathfare=' ) ;
disp( pathfar ) ;
myplot( path , Coord , pathfar ) ;