详解模拟退火算法Matlab实现
模拟退火算法概述
算法步骤
算法特点
模拟退火算法MATLAB实现
【例1】一元-多元函数优化
【例2】TSP问题
模拟退火算法概述
模拟退火算法(Simulated Annealing，简称SA)的思想最早是由Metropolis等提出的。其出发点是基于物理中固体物质的退火过程与一般的组合优化问题之间的相似性。模拟退火法是一种通用的优化算法，其物理退火过程由以下三部分组成:
p加温过程。其目的是增强粒子的热运动，使其偏离平衡位置。当温度足够高时，固体将熔为液体，从而消除系统原先存在的非均匀状态。-p
p等温过程。对于与周围环境交换热量而温度不变的封闭系统，系统状态的自发变化总是朝自由能减少的方向进行的，当span style="color:#f33b45;"自由能达到最小时，系统达到平衡状态-span。-p
p冷却过程。使粒子热运动减弱，系统能量下降，得到晶体结构。-p
p加温过程相当于对算法span style="color:#f33b45;"设定初值-span，等温过程对应算法的span style="color:#f33b45;"Metropolis抽样过程-span，冷却过程对应span style="color:#f33b45;"控制参数的下降-span。这里能量的变化就是span style="color:#f33b45;"目标函数-span，我们要得到的span style="color:#f33b45;"最优解就是能量最低态-span。其中Metropolis准则是SA算法收敛于全局最优解的关键所在，Metropolis准则以一定的概率接受恶化解，这样就使算法跳离局部最优的陷阱。-p
pSA算法的Metropolis准则允许接受一定的恶化解，具体来讲，是span style="color:#f33b45;"以一定概率来接受非最优解-span。举个例子，相当于保留一些“潜力股”，使解空间里有更多的可能性。对比轮盘赌法，从概率论来讲，它是对非最优解给予概率0，即全部抛弃。-p
p模拟退火本身是求一个最小值问题，但可以转化为求最大值问题，只需要对目标函数加个负号或者取倒数。-p
算法步骤
初始化:取初始温度T0足够大，令T = T0，任取初始解S1。
对比GA、ACA、PSO之类的群优化算法，需要在解空间中产生多个个体，再继续寻优。而SA算法只需要一个点即可。
对当前温度T，重复第(3)~(6)步。
对当前解S1随机扰动产生一个新解S2。
此处随机的扰动没有定义。结合实际例子做选择。
计算S2的增量df = f(S2) - f(S1)，其中f(S1)为S1的代价函数。
代价函数相当于之前群优化算法中讲的适应度函数。
若df 0，则接受S2作为新的当前解，即S1 = S2；否则，计算S2的接受概率exp(-df-T)， T是温度。随机产生(0,1)区间上均匀分布的随机数rand，若exp(-df-T) rand，也接受S2作为新的当前解S1 = S2，否则保留当前解S1。
这是SA算法的核心部分，即有一定概率接受非最优解。
如果满足终止条件Stop，则输出当前解S1为最优解，结束程序，终止条件Stop通常取为在连续若干个Metropolis链中新解S2都没有被接受时终止算法或者是设定结束温度。否则按衰减函数衰减T后返回第(2)步，即被接受的新的解一直在产生，则我们要对问题进行降温，使得非最优解被接受的可能不断降低，结果愈发收敛于最优解。
算法特点
与遗传算法、粒子群优化算法和蚁群算法等不同，模拟退火算法不属于群优化算法，不需要初始化种群操作。
收敛速度较慢。因为1）它初始温度一般设定得很高，而终止温度设定得低，这样才符合物体规律，认为物质处于最低能量平衡点；2）它接受恶化解，并不是全程都在收敛的过程中。这一点可以类比GA中的变异，使得它不是持续在收敛的，所以耗时更多一些。
温度管理（起始、终止温度）、退火速度（衰减函数）等对寻优
结果均有影响。比如T的衰减速度如果太快，就会导致可能寻找不到全局最优解。
模拟退火算法MATLAB实现
MATLAB自带模拟退火算法工具箱。在本文中，我们采用自带工具箱来解决函数优化的问题，再使用自己编写的程序来解决一个TSP问题。
要使用自带的工具箱，我们先安装：
【例1】一元-多元函数优化
我们要找：
一元函数：x = [1,2]范围内 y = sin(10*pi*x) - x 的极值
二元函数：在x,y都是[-5,5]范围内找z = x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20 的极值
上面是我们提前获得的ground-truth（用main.m，代码在下文）。下面我们假装不知道这个结果的，用模拟退火方法来搜索：首先定义我们在SA算法中需要用到的代价函数fitness.m：
function fitnessVal = fitness( x )
%一元函数优化：
fitnessVal = sin(10*pi*x) - x; %求最小值
% fitnessVal = -1 * sin(10*pi*x) - x; 用模拟退火求最大值，可以加个负号或者弄个倒数！
%二元函数优化：
% fitnessVal = -1 * (x(1)^2 + x(2).^2 - 10*cos(2*pi*x(1))
- 10*cos(2*pi*x(2)) + 20);
主程序main.m，用于我们直观看图先获得ground-truth：
%% I. 清空环境变量
clear all
%% II. 一元函数优化
x = 1:0.01:2;
y = sin(10*pi*x) .- x;
plot(x,y,'linewidth',1.5)
ylim([-1.5, 1.5])
xlabel('x')
ylabel('y')
title('y = sin(10*pi*x) - x')
% 1. 标记出最大值点
[maxVal,maxIndex] = max(y);
plot(x(maxIndex), maxVal, 'r*','linewidth',2)
text(x(maxIndex), maxVal, {[' X: ' num2str(x(maxIndex))];[' Y: ' num2str(maxVal)]}) % 2. 标记出最小值点
[minVal,minIndex] = min(y);
plot(x(minIndex), minVal, 'ks','linewidth',2)
text(x(minIndex), minVal, {[' X: ' num2str(x(minIndex))];[' Y: ' num2str(minVal)]})
