近五年出现的建模分析算法（模拟退火算法与数学建模）
近五年出现的建模分析算法（模拟退火算法与数学建模）在温度T，分子停留在状态r满足Boltzmann概率分布：4、数学表述热力学中的退火现象指物体逐渐降温时发生的物理現象：温度越低，物体的能量状态越低，到达足够的低点时，液体开始冷凝与结晶，在结晶状态时，系统的能量状态最低。缓慢降温（退火，annealing）时，可达到最低能量状态；但如果快速降温（淬火，quenching），会导致不是最低能态的非晶形。大自然知道慢工出细活：缓缓降温，使得物体分子在每一温度时，能够有足够时间找到安顿位置，则逐渐地，到最后可得到最低能态，系统最稳定。模仿自然界退火現象而得，利用了物理中固体物质的退火过程与一般优化问题的相似性，从某一初始温度开始，伴随温度的不断下降，结合概率突跳特性在解空间中随机寻找全局最优解。
一、前言模拟退火(Simulated Annealing，SA)算法是局部搜索算法的扩展。它源于对固体退火过程的模拟，采用Metropolis接受准则，并用一组称为冷却进度表的参数控制算法进程，使算法在多项式时间里给出一个近似最优解．
SA算法也是一种通用的随机搜索算法，它可用于解决众多的组合优化问题。当待解决的问题复杂性较高且规模较大时，在对问题的领域知识所知甚少的情况下，采用SA算法最合适。因为SA 算法不像其它确定性的启发式算法那样，需要依赖问题的领域知识来提高算法的性能。但从另一方面来说，如果已知有关待解问题的一些知识后，SA算法却无法充分利用它们，这时SA算法的优点就成了缺点。如何把传统的启发式搜索方法与SA这样的随机捜索算法结合起来，这是一个有待研究的十分有意义的课题。
SA算法在求解规模较大的实际问题时，往往存在着收敛速率慢的缺点。为此人们对SA算 法提出了各种各样的改进，它们包括并行SA算法、快速SA算法和对SA算法中各个函数和参数的重新设计等。
二、模拟退火算法及模型1、算法的提出模拟退火算法最早的思想由Metropolis等（1953）提出，1983年Kirkpatrick等将其应用于组合优化。
2、算法的目的- 解决NP复杂性问题；
- 克服优化过程陷入局部极小；
- 克服初值依赖性。
热力学中的退火现象指物体逐渐降温时发生的物理現象：温度越低，物体的能量状态越低，到达足够的低点时，液体开始冷凝与结晶，在结晶状态时，系统的能量状态最低。缓慢降温（退火，annealing）时，可达到最低能量状态；但如果快速降温（淬火，quenching），会导致不是最低能态的非晶形。
大自然知道慢工出细活：缓缓降温，使得物体分子在每一温度时，能够有足够时间找到安顿位置，则逐渐地，到最后可得到最低能态，系统最稳定。
模仿自然界退火現象而得，利用了物理中固体物质的退火过程与一般优化问题的相似性，从某一初始温度开始，伴随温度的不断下降，结合概率突跳特性在解空间中随机寻找全局最优解。
4、数学表述
在温度T，分子停留在状态r满足Boltzmann概率分布：
在同一个温度T，选定两个能量E1<E2，有
Boltzman概率分布告诉我们：
（1）在同一个温度，分子停留在能量小状态的概率大于停留在能量大状态的概率。
（2）温度越高，不同能量状态对应的概率相差越小；温度足够高时，各状态对应概率基本相同。
（3）随着温度的下降，能量最低状态对应概率越来越大；温度趋于0时，其状态趋于1。
模拟退火算法基本思想：在一定温度下，搜索从一个状态随机地变化到另一个状态；随着温度的不断下降直到最低温度，搜索过程以概率1停留在最优解。
模拟退火算法对应了一个马尔可夫链；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率1收敛于全局最优解的算法。
5、Metropolis准则（1953）- 以概率接受新状态固体在恒定温度下达到热平衡的过程可以用Monte Carlo方法（计算机随机模拟方法）加以模拟，虽然该方法简单，但必须大量采样才能得到比较精确的结果，计算量很大。
Metropolis准则在温度T下，当前状态i → 新状态j若Ej<Ei，则接受 j 为当前状态；否则，
若概率 p=exp[-(Ej-Ei)/kT] 大于[0 1)区间的随机数，则仍接受状态 j 为当前状态；
若不成立则保留状态 i 为当前状态。
在高温T下，可接受与当前状态能量差△E=Ej-Ei较大的新状态，类比组合优化问题，增大跳出局部最优解的机会；
在低温T下，只接受与当前状态能量差较小△E=Ej-Ei的新状态，类比组合优化问题，保证加快收敛到全局最优解。
6、组合优化与金属退火相似性比较|
组合优化问题 |
金属物体退火 |
|
解 |
粒子状态 |
|
最优解 |
能量最低的状态 |
|
设定初温 |
熔解过程 |
|
Metropolis抽样过程 |
等温过程 |
|
控制参数的下降 |
冷却 |
|
目标函数 |
能量 |
理论上是用一个马尔可夫链描述模拟退火算法的变化过程，因此具有全局最优性。实际应用中的模拟退火算法是一个启发式算法，它有诸多的参数需要调整，如起始温度，温度下降的方案，固定温度时的迭代长度及终止规则等，这样需要人为地调整。人为的因素，如对问题的了解，参数和规则的搭配等，会造成计算结果的差异。解决这个矛盾的方法主要通过大量的数值模拟计算，从中选择比较好的参数搭配。
1、控制参数初始温度 t0 的选取起始温度 t0 应保证平稳分布中每一状态的概率相等，应让初始接受率接近1。
初始温度t0数值计算算法：
实验表明，初始温度越高，得到高质量解的概率越大，但算法计算的时间越长。兼顾算法求解的质量和执行效率，选择初温的方法常有:随机产生一组状态，确定最大值和最小值间的差值d，利用函数得到初温，如 T0= dp，其中p为初始接受概率; 或均匀抽样一组状态，以各状态的目标值的方差为初温。
2、状态产生函数状态产生函数要解决如下问题：当前解Xi的邻域N(Xi)怎么选择？新解（新状态）Xj怎么从当前解Xi(当前状态)产生？
答：没有固定的准则。不同的问题产生新解的方法灵活多样，如TSP问题、书店买书问题等，但要遵循的基本原则是：新解的产生要用到原来搜索过程的信息，即新解与旧解之间存在关联，这是启发式算法的核心思想。
这里我们介绍几种求新解的方法:
解n元函数最大/最小值问题产生新解的方法MATLAB中内置函数的方法:
对解空间按城市编号进行编码，设初始解为s = {2，8，3，5，9，1，4，6，7}。
当前状态s的状态产生函数N(s) 即邻域可以定义为以下三种的任意一种方式：
（1）互换操作，随机交换两个城市的顺序；
初始解s的邻域N(s) = { {8，2，3，5，9，1，4，6，7} {3，8，2，5，9，1，4，6，7}，...，{2，8，1，5，9，3，4，6，7}，...，{8，2，3，5，9，1，4，7，6} }
其中一个新解s’ = {2，8，1，5，9，3，4，6，7}。 如图4-1所示。
图4-1
（2）逆序操作，两个随机位置间的城市逆序；
同理可得邻域N(s) 其中一个新解s’ = {2，8，3，4，1，9，5，6，7}。如图4-2所示。
图4-2
（3）插入操作，随机选择某点插入某随机位置。
同理可得邻域N(s) 其中一个新解s’ = {2，3，5，9，9，1，4，6，7}。如图4-3所示。
图4-3
bool型问题求新解对解空间进行编码，当前解为s = 1001，当将邻域动作定义为翻转其中一个bit时，得到当前解的邻域N(s)={0001 1101 1011 1000}，可任选一新解s’=0001。同理，当将邻域动作定义为互换相邻bit时，得到当前解的邻域P = N(s)={0101 1001 1010} 可任选一新解s’=1001。
4皇后问题求新解如图4-4所示：
图4-4
对解空间进行编码，当前解s=(2 4 1 3)，表示第1行第2列、第2行第4列、第3行第1列、第4行第3列各有一个皇后。邻域动作定义为棋盘上任意两个皇后的所在行或列进行交换，即s中任意两个元素交换位置。其邻域N(s)= {(4 2 1 3) (1 4 2 3) (3 4 1 2) (2 1 4 3) (2 3 1 4) (2 4 3 1)}。 可任选一新解s’=(4 2 1 3)。
3、接受函数---Metropolis准则Aij(t)称为接受函数或接受概率， Aij(t)表示退火温度为t时状态 i 产生 j 后，接受 j 的概率。
高温下可接受与当前状态能差较大的新状态为重要状态，而在低温下只能接受与当前状态能差较小的新状态为重要状态，这与不同温度下热运动的影响完全一致。在温度趋于零时，就不能接受任一 f(j) >f(i)的新状态 j了。
上述接受新状态的准则称为Metropolis准则，相应的算法称为Metropolis算法，这种算法的计算量显著减少。
4、温度衰减函数1> 最简单的控制参数衰减函数（指数降温）：
其中 0< α <1，是一个接近1的常数。α越接近1温度下降得越慢．这个衰减函数被许多研究者采用，他们选用的α值是0.5－0.99。这种方法简单易行，很受应用人员的欢迎。它的每一步以相同的比率下降。
2> 下降长度相等函数：
其中t0为初始温度，L为算法温度下降的总次数。这一下降方法的优点是易于操作，而且可以简单的控制温度下降的总步数。它的每一步下降长度相等，只适用于以迭代次数为停止准则的冷却进度表。
3> 经典的模拟退火衰减函数（对数降温）为：
4> 可以选择其它的衰减函数：
衰减函数的选择要权衡考虑。由于退火必须“徐徐”降温，T也必须缓慢衰减，才能确保模拟退火算法最终趋于组合优化问题的整体最优解集。但还要考虑避免过长的马尔可夫链和迭代次数增加。
5、内循环终止准则内循环终止准则，就是每一温度的迭代长度规则，即马氏链长度 Lk 的选取。
Lk值给定算法进程在第k个马氏链中进行的变换数，有限序列Lk则规定了算法进程搜索的解范围．由于多数组合优化问题的解空间规模随问题规模n呈指数型增大，为使算法最终解的质量有所保证，理应建立 Lk与n间的某种关系。指数型的关系显然是不切合实际的，因而Lk通常取为问题规模n 的一个多项式函数。
1> 固定长度
在每一温度，迭代相同的步数，步数的选取同问题的大小有关，通常采用与邻域大小直接相关的规则．如TSP中的2-opt，邻域大小是n*(n-1)/2。(n是城市数)，在同一温度，可取Lk=n n^2 100n n*(n-1)/2。
2> 由接受和拒绝的比率来控制迭代步数
当温度很高时，每一个状态被接受的频率基本相同，而且几乎所有的状态都被接受。此时，在同一温度应使迭代的步数尽量小。当温度渐渐变低时，越来越多的状态被拒绝。如果在此温度的迭代太少，则可能造成过早地陷入局部最优状态，比较直观和有效的方法是随着温度的下降，将同一温度的迭代步长增加．实现的一种方法是给定一个充分大的步长上限U和一个接受次数指标R，当接受次数等于R时，在此温度不再迭代而使温度下降，否则，一直迭代到上限步数．实现的第二种方法是给定一个接受比率指标R，迭代步长上限U和下限 L，每一温度至少迭代 L步且记录同一温度迭代的总次数和被接受的次数，当迭代超过 L步时，若接受次数同总次数的比率不小于R时，在这一温度不再迭代而开始温度下降，否则，一直迭代到上限步数。同样可以用拒绝次数为指标类似上面的讨论得到一些控制迭代步数的规则。
3> 概率控制法
6、外循环终止准则模拟退火算法从初始温度开始，通过在每一温度的迭代和温度的下降，最后达到终止原则而停止。尽管有些原则有一定理论的指导，终止原则大多是直观的。
1> 零度法（即tf 充分小）模拟退火算法的最终温度为零。因而最为简单的原则是：给定一个比较小的正数ε ，当温度tk ≤ε时，算法停止，表示已经达到最低温度。
2> 循环总数控制法
总的温度下降次数为一定值L，当温度迭代次数达到L时，停止运算。这一原则可分为两类，一类是整个算法的总迭代步数为一固定数。它表示各温度时马氏链迭代数（内循环）的总和为一个给定的数。这样很容易计算算法的复杂性，但需要合理分配内循环的长度和温度下降的次数。
3> 基于不改进规则的控制法以算法进程所得到的某些近似解为衡量标准，判断算法当前解的质量是否持续得到明显提高，从而确定是否终止算法。Kirkpatrick等人选取的停止准则是，在若干个（取s个如s=1 s=2等）相继的马氏链中解无任何变化（含优化或恶化）就终止算法。
或在一个温度和给定的迭代次数内没有改进当前的局部最优解，则停止算法。模拟退火的一个基本思想是跳出局部最优解。直观的结论是在较高的温度没能跳出局部最优解，则在低的温度跳出局部最优解的可能也比较小。由此产生上面的停止原则。
4> 接受概率控制法
给定一个指标 Xf > 0是一个比较小的数，除当前局部最优解以外，其他状态的接受概率都小于Xf 时，停止运算。记录当前局部最优解，给定一个固定的迭代次数，当在规定的次数里没有离开局部最优解或每一次计算的接受概率都小于Xf ，则在这个温度停止计算。
5> 邻域法
若采用产生概率
和接受概率
且设 f0 和 f1 分别为一个邻域内的局部最优和次最优值，当满足
时，其中N为邻域的大小，从局部最优到次优的接受概率满足(*)，而从局部最优到其他费用更高的状态的接受概率更小。直观的想法是邻域中每次至少有一个状态被接受，但当满足(*)时，除局部最优解以外状态的接受概率都小于邻域总点平均数，此时可以认为从局部最优解转移到其他状态的可能性很小，因此停止。通过(*)可得终止温度
6> 终止温度的精细估计
五、冷却进度表定义一个冷却进度表应当规定下述参数：
- 控制参数t的初值t0 ，即初始温度的选取。
- 控制参数t的衰减函数，即温度下降的规则。
- 马氏链的长度 Lk，即每一温度马氏链的迭代长度。
- 控制参数t的终值tf，即停止准则。
模拟退火算法依据Metropolis准则接受新解，因此除接受优化解外，还在一个限定范围内接受恶化解，这正是模拟退火算法与其它局部搜索算法的本质区别所在。开始时 t 值大，可能接受较差的恶化解；随着 t 的减小，只能接受较好的恶化解；最后在 t 值趋于零时，就不再接受任何恶化解了。这就使算法既可以从局部最优的陷阱中跳出，更有可能求得组合优化问题的整体最优解；又不失简单性和通用性。因此对大多数组合优化问题而言，模拟退火算法要优于其它的局部搜索算法。
七、案例分析模拟退火算法的解题思路是从一给定解开始，使用产生器和接受准则，不断把目前结构的解转为临近结构的解，持续进行“产生新解—判断—接受或舍弃”的迭代过程，有三个重要参数分别为：生成函数、接受函数（与当前解和新解的目标函数的差值有关）和温度更新函数，使算法在多项式时间里给出一个近似最优解。
例子1 旅行商问题(TSP Traveling Salesman Problem)一个商人欲到n=30个城市推销商品，每两个城市i和j之间的距离为dij，如何选择一条道路使得商人每个城市走一遍后回到起点且所走路径最短？
如图7-1所示
|
城市编号 |
城市坐标 |
城市编号 |
城市坐标 |
|
1 |
41 94 |
16 |
25 38 |
|
2 |
37 84 |
17 |
24 42 |
|
3 |
54 67 |
18 |
58 69 |
|
4 |
25 62 |
19 |
71 71 |
|
5 |
7 64 |
20 |
74 78 |
|
6 |
2 99 |
21 |
87 76 |
|
7 |
68 58 |
22 |
18 40 |
|
8 |
71 44 |
23 |
13 40 |
|
9 |
54 62 |
24 |
82 7 |
|
10 |
83 69 |
25 |
62 32 |
|
11 |
64 60 |
26 |
58 35 |
|
12 |
18 54 |
27 |
45 21 |
|
13 |
22 60 |
28 |
41 26 |
|
14 |
83 46 |
29 |
44 35 |
|
15 |
91 38 |
30 |
4 50 |
图7-1
解：冷却进度表参数：
初始温度T0 = 280
温度衰减系数α = 0.92
内循环的迭代次数（每个温度下的马尔可夫链长）Lk = 100n (n=30)
外循环的迭代次数或停止准则：在s个相继的马氏链中解无任何变化就终止算法．（如可取s=1 即前后两次外循环解无任何变化）。
目标函数，即为访问所有城市的路径长度或称为代价函数，须求其最小值，选为
新解的产生：任选访问的序号 u 和 v ，交换u 和 v的位置，即可产生新路径为（设 u < v ）
代价（目标）函数差:
第1次外循环：
图7-2
第5次外循环：
图7-3
第50次外循环：
图7-4
第100次外循环：
图7-5
最终搜索结果：
图7-6
例子2 求2个变量函数的最大值求函数y = 21.5 x1 * sin(4π * x1) x2 * sin(20π * x2)的最大值，其中x1 ∈ [-3 12.1]， x2 ∈ [4.1 5.8]。
[注：这个函数含有非常多的局部极大值，怎么设置合适的退火参数需要不断尝试。]
解：冷却进度表参数：
初始温度T0 = 100
温度衰减系数α = 0.95
内循环的迭代次数（每个温度下的马尔可夫链长）Lk = 400
外循环的迭代次数300 即终止值 tf= (α^300)*T0=(0.95^300)*100
x的参数分别为x1、x2 参数个数为2
x的下界x_lb(2) = {-3 4.1};
x的上界 x_ub(2) = {12.1 5.8};
rand(1)表示产生[0 1)的随机数
初始化解: x0(i) = x_lb(i) (x_ub(i) - x_lb(i)) * rand(1) i=1 2
目标函数题目已给，目标函数差Δf 容易得到。
生成新解：
生成[0 1)的随机数: y(i) = rand(1) i=1 2
根据新解的产生规则计算: z(i) = y(i) / sqrt(∑(y(i) * y(i)) ) i=1 2
根据新解的产生规则计算: x_new(i) = x(i) z(i)*T i=1 2
新解边界调整：
超出下界调整x_new(i) = rand(1)*x_lb(i) (1- rand(1))*x(i) i=1 2
超出上界调整x_new(i) = rand(1)*x_ub(i) (1-rand(1))*x(i) i=1 2
两种计算机实现方式：
1> 用Python语言实现动画与模拟退火过程
略
2> 基于MATLAB内置优化工具箱的实现
（1）在MATLAB主界面中选择“APP” - “Optimization”，进入优化工具箱；
（2）在左侧面板的“Solver”中选择“simulannealbnd - Simulated annealing algorithm”，即模拟退火算法；
如下MATLAB代码仅供参考：
%% SA 模拟退火: 以求二元函数的最大值为例
% 求解函数y = 21.5 x1*sin(4*pi*x1) x2*sin(20*pi*x2)
% 在x1∈[-3 12.1] x2∈[4.1 5.8]内的最大值（可绘图动画演示搜索过程）
clear; clc
%% 绘制函数的图形
x1 = -3: 0.1: 12.1;
x2 = 4.1: 0.1: 5.8;
[x1 x2] = meshgrid(x1 x2);
y = 21.5 x1.*sin(4*pi*x1) x2.*sin(20*pi*x2); % 注意这里要点乘的地方
mesh(x1 x2 y)
xlabel('x1'); ylabel('x2'); zlabel('y'); % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比，使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on % 不关闭图形，继续在上面画图
%% 参数初始化
narvs = 2; % 变量个数
T0 = 100; % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 300; % 外层循环的最大迭代次数
Lk = 400; % 每个温度下的迭代次数
alfa = 0.95; % 温度衰减系数
x_lb = [-3 4.1]; % x的下界
x_ub = [12.1 5.8]; % x的上界
%% 随机生成一个初始解x0，并将当前搜索点记为S
x0 = zeros(1 narvs);
for i = 1: narvs
x0(i) = x_lb(i) (x_ub(i)-x_lb(i))*rand(1);
end
S = x0;
h = scatter3(S(1) S(2) Obj_fun2(S) '*r'); % scatter3是绘制三维散点图的函数（这里返回h是为了得到图形的句柄，未来我们对其位置进行更新）
%% 初始化用来保存中间结果的x和y的取值
iter_x = S;
iter_y = Obj_fun2(S);
%% 模拟退火过程
for iter = 1 : maxgen % 外循环 我这里采用的是指定最大迭代次数
for i = 1 : Lk % 内循环（迭代次数即马尔科夫链长度），在每个温度下开始迭代
% 计算当前解的函数值
y_now = 21.5 S(1)*sin(4*pi*S(1)) S(2)*sin(20*pi*S(2));
% 随机产生一个新解
y = randn(1 narvs); % 生成1行narvs列的N(0 1)随机数
z = y / sqrt(sum(y^2)); % 根据新解的产生规则计算z
x_new = S z*T; % 根据新解的产生规则计算x_new的值
% 判断新解的位置是不是超出了定义域，如果超出了定义域，就对其进行调整
for j = 1: narvs
if x_new(j) < x_lb(j)
r = rand(1);
x_new(j) = r*x_lb(j) (1-r)*S(j);
elseif x_new(j) > x_ub(j)
r = rand(1);
x_new(j) = r*x_ub(j) (1-r)*S(j);
end
end
% 计算新解的函数值
y_new = 21.5 x_new(1)*sin(4*pi*x_new(1)) x_new(2)*sin(20*pi*x_new(2));
% 更新搜索点的位置
if y_new > y_now % 如果新解函数值大于当前解的函数值
S = x_new; % 搜索点移动到新解的位置
iter_x = [iter_x; x_new]; % 将新找到的x1添加到中间结果中
iter_y = [iter_y; y_new]; % 将新找到的y1添加到中间结果中
else
p = exp(-abs(y_new - y_now)/T); % 根据Metropolis准则计算一个概率
if rand(1) < p % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
S = x_new; % 搜索点移动到新解的位置
end
end
% 内循环结束！注意内循环不需要对温度进行衰减
end
T = alfa * T; % 温度下降
pause(0.05) % 暂停一段时间后(单位：秒)再接着画图
h.XData = S(1); % 更新散点图句柄的x轴的数据（此时解的位置在图上发生了变化）
h.YData = S(2); % 更新散点图句柄的y轴的数据（此时解的位置在图上发生了变化）
h.ZData = Obj_fun2(S); % 更新散点图句柄的z轴的数据（此时粒子的位置在图上发生了变化）
end
%% 查找搜索全过程的最大值
[best_y ind] = max(iter_y); % 找到最大的y的值，以及其在向量中的下标
best_x = iter_x(ind :); % 根据下标找到此时x的值
disp('最佳的位置是：'); disp(best_x)
disp('此时最优值是：'); disp(best_y)
%% 在图上标注出搜索到的最大值点
h.XData = []; h.YData = []; h.ZData = [] % 将原来的散点删除
scatter3(best_x(1) best_x(2) best_y '*r'); % 在最大值处重新标上散点
title(['模拟退火找到的最小值为' num2str(best_y)]) % 加上图的标题例子3
[注：近似最优解x=[-0.0282 0.0046 -0.0158 0.0265 0.0345 0.0436 -0.0467 0.0006 0.0179 -0.0282]，函数f(x)的最小值为0.008156。]
解：冷却进度表参数：
初始温度T0 = 100
温度衰减系数α = 0.998
内循环的迭代次数（每个温度下的马尔可夫链长）Lk = 200
外循环的迭代次数：容差1e-8与终止值tf= 0.001
x的参数个数为10
x的下界x_lb = -20;
x的上界 x_ub = 20;
rand(1)表示产生[0 1)的随机数
初始化解: x0(i) = x_lb (x_ub - x_lb) * rand(1) i=1 2 ... 10
目标函数题目已给，目标函数差Δf 容易得到。
生成新解：
x_new[i]=x[i] S*(x_lb (x_ub - x_lb) * rand(1)) i=1 ... 10 步长因子S=0.01。
如下MATLAB代码仅供参考：
%初始化
clear all; %清除所有变量
close all; %清图
clc; %清屏
D=10; %变量维数
Xs=20; %上限
Xx=-20; %下限
%冷却表参数
L=200; %马尔科夫链长度
K=0.998; %衰减参数
S=0.01; %步长因子
T=100; %初始温度
YZ=1e-8; %容差
P=0; %Metropolis过程中总接受点
%随机选定初值设定
PreX=rand(D 1)*(Xs-Xx) Xx;
PreBestX=PreX;
PreX=rand(D 1)*(Xs-Xx) Xx;
BestX=PreX;
%每迭代一次退火一次（降温），直到满足迭代条件为止
deta=abs(func1(BestX)-func1(PreBestX));
while (deta>YZ)&&(T>0.001)
T=K*T;
%在当前温度T下迭代次数
for i=1:L
%在此点附近随机选下一个点
NextX=PreX S*(rand(D 1)*(Xs-Xx) Xx);
%边界条件处理
for ii=1:D
if NextX(ii)>Xs|NextX(ii)<Xx
NextX(ii)=PreX(ii) S*(rand*(Xs-Xx) Xx);
end
end
%是否全局最优解
if(func1(BestX)>func1(NextX))
%保留上一个最优解
PreBestX=BestX;
%此为新的最优解
BestX=NextX;
end
%Metropolis过程
if(func1(PreX)-func1(NextX)>0)
%接受新解
PreX=NextX;
P=P 1;
else
changer=-1*(func1(NextX)-func1(PreX))/T;
p1=exp(changer);
%接受较差的解
if p1>rand
PreX=NextX;
P=P 1;
end
end
trace(P 1)=func1(BestX);
end
deta=abs(func1(BestX)-func1(PreBestX));
end
disp('最小值在点：');
BestX
disp('最小值为：');
func1(BestX)
figure
plot(trace(2:end))
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')
%适应度函数
function result=func1(x);
summ=sum(x.^2);
result=summ;例子4
解：这是一个非线性量的非线性目标函数且带有非线性约束的优化问题，用理论方法求解有一定的困难，用模拟退火算法则容易得多。
实际求解时，最简单的方法是每次随机产生变量 x、y、z 的一组值，当满足约束时再求目标函数的值，并按算法重复迭代，直至算法终止．但对本例这种带有等式约束的问题，在相当多的重复试验中也很难产生出满足约束的变量值组。为此可先将约束中的等式处理为
y 的取值范围，由约束条件可得
因此可在此范围内随机产生y的值，然后对(*)式求出的 z 和x判断是否非负即可．
冷却进度表参数：
初始温度T0 = 0.1
温度衰减系数α = 0.5
内循环的迭代次数（每个温度下的马尔可夫链长）Lk = 40
外循环的迭代次数或停止准则：在s个相继的马氏链中解无任何变化就终止算法．（如可取s=1 即前后两次外循环解无任何变化）。
初始解选为 x = y = z = 0 f =0（最终的近似最优解与初始解的选择无关）。
新解x，y，z产生规则：在y的约束范围随机取值y，代入(*) 可求得x z。
目标函数题目已给，目标函数差Δf 容易得到。
求出的近似最优解为：当x=21.5905 y =0.562582 z =5.77765，fmax=345.2348。
理论最优解为：当x=21.5912，y=0.567274 z=5.77652 时，fmax=345.2350。
例子5 0-1背包问题（ZKP）给定一个可装重量M的背包n件物品，物品i的重量和价值分别为和 i=1 n。要选若干件物品装人背包，使其价值之和最大。
解：模拟退火算法描述如下：
① 解空间ZKP是一个有约束的优化问题。为此，我们限定解空间为所有可行解的集合，即
其中xi = 1表示物品i被选择装入背包。初始解选为(0 … 0)。
② 目标函数求最大值的代价函数
但必须满足约束
随机选取物品i。若i不在背包中，则将其直接装入背包，或同时从背包中随机取出另一件物品j。若i已在背包中，则将其取出，同时随机装入另一件物品j。即
根据产生新解的三种可能，伴随的背包代价差为
为判定解的可行性，还需要求出对应的背包重量差为
其中∆m为当前背包重量m的增重
⑤ 接受准则是扩充了可行性测定的马尔可夫链准则
给定带权图G={V E} 其中V={v1 ... vn}为顶点集合，E为边集，加权矩阵为W=[w(i j)]。要将V划分为子集V0和V1，使E中所有顶点分属V0和V1的权之和最大。
① 解空间可表为V的所有划分为两子集和V0分划V1的集合，即
S={δ|δ为将V划分为V0和V1的一个和分划}
其中分划 δ(i) = j i=1 ... n; j=0 1
表示顶点vi ∈ Vj 。初始解选为δ(i) = 0 i=1 ... n。
② 目标函数直接取为分划δ所得得到的截量的相反数
需求其最小值。即相应截量的最大值。
③ 新解的产生任选顶点u∈V，将其从目前所在的子集移到另一个子集中去，即:
δ(u) = 1-δ(u)
④ 目标函数差
根据目标函数和产生新解的方法可知，相应于目标函数差为
例如：6个顶点的带权图
第1步，初始化：
V0={v1 v2 v3 v4 v5 v6} δ(i)=0 i=1 ... 6
f(δ) =0
第2步，任选v5 将v5从V0转移到V1 V0={v1 v2 v3 v4 v6} V1={v5}。 δ(5) =1。
f(δ) = -(w(5 6) w(5 1) w(5 7) w(5 4)) = -45
∆f=-45-0 =-45
第3步，任选v4 将v4从V0转移到V1 V0={v1 v2 v3 v6} V1={v5 v4}。 δ(5) =1 δ(4) =1。
f(δ) = -(w(5 6) w(5 1) w(5 7) w(4 3) w(4 2) w(4 7)) = -51。
∆f=-51-(-45) =-6
或者∆f=w(4 5) - (w(4 7) w(4 2) w(4 3)) = 13-19=-6
冷却进度表参数：
初始温度T0 = 0.8 (参考范围0.5-2)
温度衰减系数α = 0.95 （参考范围 大于0.9)
内循环的迭代次数（每个温度下的马尔可夫链长）Lk = 1000 (参考值 500，1000等）
外循环的迭代次数或停止准则：在s个相继的马氏链中解无任何变化就终止算法．（如可取s=1 即前后两次外循环解无任何变化）。
例子7 独立集问题（ISP）
模拟退火算法的求解形式为：
解空间
求解的模拟退火算法为
解空间
为了用模拟退火算法求解，首先注意到，若图G的顶点个数为n，则图G的最小着色的上界为n。此外，GCP也是一个有约束的优化问题，我们同样使用惩函数将其转化为无约束问题。
解空间
(1) 目前人们已发现千余个NP完全问题，这些NP完全问题中有许多可以用模拟退火算法求解，以上给出的仅仅是其中很小的一个子集。
(2) 模拟退火算法是一种随机性的近似算法，其效率和所得解的质量依赖于冷却进度表的选取，还受到所用的随机数序列及表达问题的数学模型的影响。在实际应用时应根据问题的具体情况适当选择，以尽量提高算法的效率和解的质量。
(3) 对于某些具体问题或实例而言，模拟退火算法的性能(所需时间和所得解的质量)不一定优于目前已有的启发式算法如：禁忌算法、遗传算法、人工蚁算法和神经网络算法。但由于该算法具有灵活性、通用性和健壮性，其优越性仍是很明显的。
(4) 模拟退火算法应用非常广泛，下面列举一些研究文章。
[1] 模拟退火算法及其应用 计算机研究与发展，1990第7期，姚新，陈国良
[2] 模拟退火法及其收敛性 光电工程 第20卷 第3期，1993年6月 陈小刚，林大键，孙国良
[3]《 现代优化计算方法 》— 邢文训，谢金星
[4]《 非数值并行算法 第一册 — 模拟退火算法》— 康立山，谢云等
