模拟退火算法来源于固体退火原理，是一种基于概率的算法，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小，即在局部最优解能概率性地跳出并最终趋于全局最优。参考了金属冶炼的退火过程。

---

模拟退火的流程

假设一个人在一群山峰中的某一个位置，他想要找一个最低点，只需要一直往比自己位置低的方向走，直到周围的地势都比自己所处的位置高就行了，这就是梯度下降法，如果想找最高点反着来就行了。但是这样面对的问题就是不一定能找到整片山峰中的最低点，如果他运气好位于最低的一座山上就没问题，但他大概率不会处于这个位置，只能找到某一座山的最低点，即局部最优解

而模拟退火算法就是有一定的概率跳到另外一座山上，而且是在新解的值比当前解的值差的情况下，如果新解比当前解情况好可以当成普通梯度下降法

具体过程：

1. 随机生成一个初始值，并产生一个初始位移，求出系统因此而产生的能量变化 $\Delta E_k$
2. 若 $\Delta E_k\le 0$ 则位移生效并继续产生位移再判断；若 $\Delta E_k>0$ 则以值为 $P_k$ 的概率来判断是否生效
   $P_k$ 公式，式中 T 为温度：$$P_k=\frac{1}{1+e^{-\Delta E_k/T}}$$
3. 转第 1 步继续执行，直到达到平衡状态为止

流程图：

---

案例解读

求目标函数在约束下的最小值
$$\begin{gathered} f(x)=x_1^2+x_2^2+8 \\ 约束：x_1^2-x_2>0-x_1-x_2^2+2=0 \end{gathered}$$

下面代码采用 MatLab，把同一个文件中的代码分几段解读

初始化

sol_new2 = 1; %初始解，即 x2 的新值
sol_new1 = 2-sol_new2^2; %x1 的新值
sol_current1 = sol_new1; %x1 当前值
sol_best1 = sol_new1; %x1 最优值
sol_current2 = sol_new2; %x2 当前值
sol_best2 = sol_new2; %x2 最优值
E_current = inf; %E 当前值初始为无穷大
E_best = inf; %E 最优值初始为无穷大rand('state', sum(clock)); %初始化随机数发生器
t = 90; %初始温度
tf = 89.9; %结束温度
a = 0.99; %温度下降比例

初始化之后就可以随机生成位移并判断，降温时新温度 $t=t\ast a$

生成新值

while t >= tf %结束条件，以温度判断for r = 1:500 %退火次数%产生随机扰动，生成新解sol_new2 = sol_new2+rand*0.2;sol_new1 = 2-sol_new2^2;%检查是否满足约束条件if sol_new1^2-sol_new2 >= 0 && -sol_new1-sol_new2^2+2==0 && sol_new1>=0 &&sol_new2>=0elsesol_new2 = rand*2;sol_new1 = 2-sol_new2^2;continue;end

首先满足默认约束条件，再考虑值的范围，要把所有约束都写上

结束条件不仅要判断退火次数，还要判断温度是否足够

退火过程

        %退火过程E_new = sol_new1^2+sol_new2^2+8; %目标函数if E_new < E_current %接受准则E_current = E_new;sol_current1 = sol_new1;sol_current2 = sol_new2;if E_new < E_best%把冷却过程中最好的解保存下来E_best = E_new;sol_best1 = sol_new1;sol_best2 = sol_new2;endelseif rand < exp(-(E_new-E_current)/t) %代价函数差E_current = E_new;sol_current1 = sol_new1;sol_current2 = sol_new2;elsesol_new1 = sol_current1;sol_new2 = sol_current2;endendplot(r,E_best, '*')hold onendt = t*a; %降温
end

上述代码就是按照模拟退火的过程进行的，下一步如果是值更优则采纳，如果更差则判断是否跳动，再继续产生新值来判断，如此重复进行，直至满足条件

展示结果

如图：

可以看到在 200 多次后，值已经非常接近最优值，每次的情况都不太一样，次数不能太小

打印最优值：

disp('最优解为：')
disp(sol_best1)
disp(sol_best2)
disp('目标表达式的最小值等于：')
disp(E_best)