Navigator
Simulated Annealing
SA是一种适合求解大规模组合优化问题的算法，是一种关于NP完全类问题的有效近似算法. SA算法是基于Monte-Carlo迭代求解策略的一种随机寻优算法，算法采用Metropolis准则，并使用一组冷却进度表的参数控制算法的进程，使得算法可以在多项式时间内给出近似最优解.
Metropolis准则
算法的模拟过程通过Monte-Carlo方法进行，采样时着重选取那些有着重要贡献的样本，Metropolis等在1953年提出了重要性采样方法，即以一定的概率接受新状态. 首先给定粒子相对位置表征的初始状态
i
i
i作为固体当前的状态，能量为
E
i
E_i
Ei，摄动微小变化到一个新状态
E
j
E_j
Ej，如果
E
j
<
E
i
E_j<E_i
Ej<Ei，那么该状态就作为重要状态；如果
E
j
>
E
i
E_j>E_i
Ej>Ei，依概率判断，固体处于
i
i
i和
j
j
j的概率比值等于玻尔兹曼常数的比值，即
r
=
exp
(
E
i
−
E
j
k
T
)
r=\exp(\frac{E_i-E_j}{kT})
r=exp(kTEi−Ej)
大量迁移后，系统处于稳定状态，概率分布趋于以下吉布斯正则分布
P
i
=
1
Z
exp
(
−
E
i
k
T
)
P_i=\frac{1}{Z}\exp\bigg(\frac{-E_i}{kT}\bigg)
Pi=Z1exp(kT−Ei)
其中
P
i
P_i
Pi为处于某微观状态或其附近的概率分布；
Z
Z
Z为常数.
SA基本过程
SA由某一个较高的初温开始，利用具有概率突跳特性的Metropolis抽样策略在解空间进行随机搜索，伴随温度不断下降的重复抽样的过程，算法步骤可以描述如下：
- 给定初温 t = t 0 t=t_0 t=t0，随机产生初始状态, s = s 0 s=s_0 s=s0,令 k = 0 k=0 k=0.
- repeat
- 产生新状态 s j = G e n e r a t e ( s ) s_j=Generate(s) sj=Generate(s)
- 可以使用不同策略产生新解，一般方法是在当前解的基础上产生新解
s j = s i + Δ s , Δ s = y ( U B − L B ) s_j=s_i+\Delta s, \Delta s=y(UB-LB) sj=si+Δs,Δs=y(UB−LB)
其中， y y y为零两侧对称分布的随机数，随机数的分布由概率密度决定， U B UB UB和 L B LB LB各位参数区间的上下界. - 如果 min { 1 , exp [ − ( C ( s j ) − C ( s ) ) / t k ] } ≥ r a n d o m [ 0 , 1 ] , s = s j \min\{1, \exp[-(C(s_j)-C(s))/t_k]\}\geq random[0, 1], s=s_j min{1,exp[−(C(sj)−C(s))/tk]}≥random[0,1],s=sj.
- 退温
- 直到算法满足终止准则
- 输出搜索结果
SA的控制参数
合理选择一组控制算法进程的参数，用以逼近模拟SA算法的渐近收敛的状态，使得算法在有限时间内返回一个近似最优解，这组控制参数一般被称为冷却进度表，包括以下参数：1.控制参数 t t t的初始温度；2.控制参数 t t t的衰减函数 t k = f ( k ) t_k=f(k) tk=f(k)；3.控制参数 t t t的终值（停止准则）4.马尔科夫链的长度 L k L_k Lk，所有迭代过程称为一个马尔科夫链.
Demo:求极小值
f
(
x
)
=
(
a
−
b
x
1
2
+
x
1
4
/
3
)
x
1
2
+
x
1
x
2
+
(
−
c
+
c
x
2
2
)
x
2
2
f(x)=(a-bx_1^2+x_1^4/3)x_1^2+x_1x_2+(-c+cx_2^2)x_2^2
f(x)=(a−bx12+x14/3)x12+x1x2+(−c+cx22)x22
其中
a
,
b
,
c
a, b, c
a,b,c为常数.
%% 使用模拟退火算法求解最小值
global fitall_sa
fitall_sa = zeros(1, 1000);
a=4; b=2.1; c=4; x0=[0.5 0.5];
% ops = optimoptions(@simulannealbnd, 'ReannealInterval', 300, 'MaxIter', 1000, 'PlotFcns', @saplotbestf);
ops= optimoptions(@simulannealbnd,'OutputFcn', @myoutSA);
[x, fval, exitflag, output]=simulannealbnd(@(x)obj_func_sa5(x, a, b, c), x0, [], [], ops);
% plot
plot(1:1000, fitall_sa(1:1000));
grid on;
Reference
最优化方法及其MATLAB实现