摘要
本文简介了模拟退火的基本思想，以于模拟时的主要参数的选择根据，然后给出一个求二维函数极值的具体问题和解法 , 并给出C# 源代码。
l 概述
在管理科学、计算机科学、分子物理学和生物学以及超大规模集成电路设计、代码设计、图像处理和电子工程等科技领域中，存在大量组合优化瓿。其中许多问题如货郎担问题、图着色问题、设备布局问题以及布线问题等，至今没有找到有效的多项式时间算法。这些问题已被证明是 NP 完全问题。
1982 年 ,KirkPatrick 将退火思想引入组合优化领域，提出一种解大规模组合优化问题的算法，对 NP 完全组合优化问题尤其有效。这源于固体的退火过程，即先将温度加到很高 , 再缓慢降温 ( 即退火 ) ，使达到能量最低点。如果急速降温 ( 即为淬火 ) 则不能达到最低点 . 。
在［ 1 ］中的解释为： Simulation Annealing is a technique which can be applied to any minimisation or learning process based on successive update steps (either random or deterministic) where the update step length is proportional to an arbitrarily set parameter which can play the role of a temperature. Then, in analogy with the annealing of metals, the temperature is made high in the early stages of the process for faster minimisation or learning, then is reduced for greater stability.
即：模拟退火算法是一种能应用到求最小值问题或基本先前的更新的学习过程（随机或决定性的）。在此过程中，每一步更新过程的长度都与相应的参数成正比，这些参数扮演着温度的角色。然后，与金属退火原理相类似，在开始阶段为了更快地最小化或学习，温度被升得很高，然后才（慢慢）降温以求稳定。
l 模拟退火算法的主要思想
就函数最小值问题来说，模拟退火的主要思想是：在搜索区间（二维平面中）随机游走（即随机选择点），再以 Metropolis 抽样准则，使随机游走逐渐收敛于局部最优解。而温度即是 Metropolis 算法中的一个重要控制参数，可以认为这个参数的大小控制了随时过程向局部或全局最优解移动的快慢。
冷却参数表、领域结构和新解产生器、接受准则和随机数产生器（即 Metropolis 算法）一起构成算法的三大支柱。
l 重点抽样与 Metroplis 算法：
Metropolis 是一种有效的重点抽样法，其算法为：系统从能量一个状态变化到另一个状态时，相应的能量从 E1 变化到 E2 ，概率为 p = exp[ - (E2- E1)/kT ] 。如果 E2 < E1 ，系统接收此状态，否则，以一个随机的概率接收此或丢弃此状态。这种经常一定次数的迭代，系统会逐渐趋于一引稳定的分布状态。
重点抽样时，新状态下如果向下则接受（局部最优），若向上（全局搜索），以一定机率接受。模拟退火方法从某个初始解出发 , 经过大量解的变换后 , 可以求得给定控制参数值时组合优化问题的相对最优解。然后减小控制参数 T 的值，重复执行 Metropolis 算法，就可以在控制参数 T 趋于零时，最终求得组合优化问题的整体最优解。控制参数的值必须缓慢衰减。
其中温度是一个 Metropolis 的重要控制参数，模拟退火可视为递减控制参数什时 Metroplis 算法的迭代。开始 T 值大，可能接受较差的恶化解，随着 T 的减小，只能接受较好的恶化解，最后在 T 趋于 0 时，就不再接受任何恶化解了。
在无限高温时，系统立即均匀分布，接受所有提出的变换。 T 的衰减越小， T 到达终点的时间越长；但可使马可夫链越小，到达准平衡分布的时间越短，
l 参数的选择：
我们称调整模拟退火法的一系列重要参数为冷却进度表。它控制参数 T 的初值及其衰减函数，对应的 MARKOV 链长度和停止条件，非常重要。
一个冷却进度表应当规定下述参数：
1 ． 控制参数 t 的初值 t0;
2 ． 控制参数 t 的衰减函数；
3 ． 马尔可夫链的长度 Lk 。（即每一次随机游走过程，要迭代多少次，才能趋于一个准平衡分布，即一个局部收敛解位置）
4 ． 结束条件的选择
有效的冷却进度表判据：
一．算法的收敛：主要取决于衰减函数和马可夫链的长度及停止准则的选择
二．算法的实验性能：最终解的质量和 CPU 的时间
参数的选择：
一）控制参数初值 T0 的选取
一般要求初始值 t0 的值要充分大，即一开始即处于高温状态，且 Metropolis 的接收率约为 1 。
二）衰减函数的选取
衰减函数用于控制温度的退火速度，一个常用的函数为： T(n + 1) = K*T(n) ，其中 K 是一个非常接近于 1 的常数。
三）马可夫链长度 L 的选取
原则是：在衰减参数 T 的衰减函数已选定的前提下， L 应选得在控制参数的每一取值上都能恢复准平衡。
四）终止条件
有很多种终止条件的选择，各种不同的条件对算法的性能和解的质量有很大影响，本文只介绍一个常用的终止条件。即上一个最优解与最新的一个最优解的之差小于某个容差，即可停止此次马尔可夫链的迭代。
使用模拟退火法求函数 f(x,y) = 5sin(xy) + x2 + y2 的最小值
解：根据题意，我们设计冷却表进度表为：
即初始温度为 100
衰减参数为 0.95
马可夫链长度为 10000
Metropolis 的步长为 0.02
结束条件为根据上一个最优解与最新的一个最优解的之差小于某个容差。
使用METROPOLIS接受准则进行模拟, 程序如下
/*
* 模拟退火法求函数f(x,y) = 5sin(xy) + x^2 + y^2的最小值
* 结束条件为两次最优解之差小于某小量
*/
using System;
namespace SimulateAnnealing
{
class Class1
{
// 要求最优值的目标函数
static double ObjectFunction( double x, double y )
{
double z = 0.0;
z = 5.0 * Math.Sin(x*y) + x*x + y*y;
return z;
}
[STAThread]
static void Main(string[] args)
{
// 搜索的最大区间
const double XMAX = 4;
const double YMAX = 4;
// 冷却表参数
int MarkovLength = 10000; // 马可夫链长度
double DecayScale = 0.95; // 衰减参数
double StepFactor = 0.02; // 步长因子
double Temperature = 100; // 初始温度
double Tolerance = 1e-8; // 容差
double PreX,NextX; // prior and next value of x
double PreY,NextY; // prior and next value of y
double PreBestX, PreBestY; // 上一个最优解
double BestX,BestY; // 最终解
double AcceptPoints = 0.0; // Metropolis过程中总接受点
Random rnd = new Random();
// 随机选点
PreX = -XMAX * rnd.NextDouble() ;
PreY = -YMAX * rnd.NextDouble();
PreBestX = BestX = PreX;
PreBestY = BestY = PreY;
// 每迭代一次退火一次(降温), 直到满足迭代条件为止
do
{
Temperature *=DecayScale;
AcceptPoints = 0.0;
// 在当前温度T下迭代loop(即MARKOV链长度)次
for (int i=0;i<MarkovLength;i++)
{
// 1) 在此点附近随机选下一点
do
{
NextX = PreX + StepFactor*XMAX*(rnd.NextDouble()-0.5);
NextY = PreY + StepFactor*YMAX*(rnd.NextDouble()-0.5);
}
while ( !(NextX >= -XMAX && NextX <= XMAX && NextY >= -YMAX && NextY <= YMAX) );
// 2) 是否全局最优解
if (ObjectFunction(BestX,BestY) > ObjectFunction(NextX,NextY))
{
// 保留上一个最优解
PreBestX =BestX;
PreBestY = BestY;
// 此为新的最优解
BestX=NextX;
BestY=NextY;
}
// 3) Metropolis过程
if( ObjectFunction(PreX,PreY) - ObjectFunction(NextX,NextY) > 0 )
{
// 接受, 此处lastPoint即下一个迭代的点以新接受的点开始
PreX=NextX;
PreY=NextY;
AcceptPoints++;
}
else
{
double change = -1 * ( ObjectFunction(NextX,NextY) - ObjectFunction(PreX,PreY) ) / Temperature ;
if( Math.Exp(change) > rnd.NextDouble() )
{
PreX=NextX;
PreY=NextY;
AcceptPoints++;
}
// 不接受, 保存原解
}
}
Console.WriteLine("{0},{1},{2},{3}",PreX, PreY, ObjectFunction ( PreX, PreY ), Temperature);
} while( Math.Abs( ObjectFunction( BestX,BestY) – ObjectFunction (PreBestX, PreBestY)) > Tolerance );
Console.WriteLine("最小值在点:{0},{1}",BestX, BestY);
Console.WriteLine( "最小值为:{0}",ObjectFunction(BestX, BestY) );
}
}
}
l 结果：
最小值在点:-1.07678129318956,1.07669421564618
最小值为:-2.26401670947686
l 后记：
原来想写一系列的文章，介绍如何用C＃解数值问题，这是因为在CSDN上很少有数值计算方面的文章，所以希望能有所补充。
一开始在网上搜索模拟退火的资料并想作为 C ＃数值计算的一个例子，找不到现成的源码。后来自己实验了很久，终于将此程序写出，不敢私藏，拿出来作用模拟退火或者用 C ＃解数值算法问题的一个入门例子。
本文尽量避免太过学术化，如数学和物理名称和公式，仓促下笔，有很多地方可能讲得不是很清楚，希望各位体谅。任何问题或批评，可 EMAIL 与我： armylau2@163.com
另，模拟退火还可以应用到其它更多更复杂的问题，如“推销员问题”等组合优化问题。本例只是求一个二维函数的最小值问题，而且其冷却表参数的选择也过于简单，只能作用一个初步的入门简介，请读者注意。
本文简介了模拟退火的基本思想，以于模拟时的主要参数的选择根据，然后给出一个求二维函数极值的具体问题和解法 , 并给出C# 源代码。
l 概述
在管理科学、计算机科学、分子物理学和生物学以及超大规模集成电路设计、代码设计、图像处理和电子工程等科技领域中，存在大量组合优化瓿。其中许多问题如货郎担问题、图着色问题、设备布局问题以及布线问题等，至今没有找到有效的多项式时间算法。这些问题已被证明是 NP 完全问题。
1982 年 ,KirkPatrick 将退火思想引入组合优化领域，提出一种解大规模组合优化问题的算法，对 NP 完全组合优化问题尤其有效。这源于固体的退火过程，即先将温度加到很高 , 再缓慢降温 ( 即退火 ) ，使达到能量最低点。如果急速降温 ( 即为淬火 ) 则不能达到最低点 . 。
在［ 1 ］中的解释为： Simulation Annealing is a technique which can be applied to any minimisation or learning process based on successive update steps (either random or deterministic) where the update step length is proportional to an arbitrarily set parameter which can play the role of a temperature. Then, in analogy with the annealing of metals, the temperature is made high in the early stages of the process for faster minimisation or learning, then is reduced for greater stability.
即：模拟退火算法是一种能应用到求最小值问题或基本先前的更新的学习过程（随机或决定性的）。在此过程中，每一步更新过程的长度都与相应的参数成正比，这些参数扮演着温度的角色。然后，与金属退火原理相类似，在开始阶段为了更快地最小化或学习，温度被升得很高，然后才（慢慢）降温以求稳定。
l 模拟退火算法的主要思想
就函数最小值问题来说，模拟退火的主要思想是：在搜索区间（二维平面中）随机游走（即随机选择点），再以 Metropolis 抽样准则，使随机游走逐渐收敛于局部最优解。而温度即是 Metropolis 算法中的一个重要控制参数，可以认为这个参数的大小控制了随时过程向局部或全局最优解移动的快慢。
冷却参数表、领域结构和新解产生器、接受准则和随机数产生器（即 Metropolis 算法）一起构成算法的三大支柱。
l 重点抽样与 Metroplis 算法：
Metropolis 是一种有效的重点抽样法，其算法为：系统从能量一个状态变化到另一个状态时，相应的能量从 E1 变化到 E2 ，概率为 p = exp[ - (E2- E1)/kT ] 。如果 E2 < E1 ，系统接收此状态，否则，以一个随机的概率接收此或丢弃此状态。这种经常一定次数的迭代，系统会逐渐趋于一引稳定的分布状态。
重点抽样时，新状态下如果向下则接受（局部最优），若向上（全局搜索），以一定机率接受。模拟退火方法从某个初始解出发 , 经过大量解的变换后 , 可以求得给定控制参数值时组合优化问题的相对最优解。然后减小控制参数 T 的值，重复执行 Metropolis 算法，就可以在控制参数 T 趋于零时，最终求得组合优化问题的整体最优解。控制参数的值必须缓慢衰减。
其中温度是一个 Metropolis 的重要控制参数，模拟退火可视为递减控制参数什时 Metroplis 算法的迭代。开始 T 值大，可能接受较差的恶化解，随着 T 的减小，只能接受较好的恶化解，最后在 T 趋于 0 时，就不再接受任何恶化解了。
在无限高温时，系统立即均匀分布，接受所有提出的变换。 T 的衰减越小， T 到达终点的时间越长；但可使马可夫链越小，到达准平衡分布的时间越短，
l 参数的选择：
我们称调整模拟退火法的一系列重要参数为冷却进度表。它控制参数 T 的初值及其衰减函数，对应的 MARKOV 链长度和停止条件，非常重要。
一个冷却进度表应当规定下述参数：
1 ． 控制参数 t 的初值 t0;
2 ． 控制参数 t 的衰减函数；
3 ． 马尔可夫链的长度 Lk 。（即每一次随机游走过程，要迭代多少次，才能趋于一个准平衡分布，即一个局部收敛解位置）
4 ． 结束条件的选择
有效的冷却进度表判据：
一．算法的收敛：主要取决于衰减函数和马可夫链的长度及停止准则的选择
二．算法的实验性能：最终解的质量和 CPU 的时间
参数的选择：
一）控制参数初值 T0 的选取
一般要求初始值 t0 的值要充分大，即一开始即处于高温状态，且 Metropolis 的接收率约为 1 。
二）衰减函数的选取
衰减函数用于控制温度的退火速度，一个常用的函数为： T(n + 1) = K*T(n) ，其中 K 是一个非常接近于 1 的常数。
三）马可夫链长度 L 的选取
原则是：在衰减参数 T 的衰减函数已选定的前提下， L 应选得在控制参数的每一取值上都能恢复准平衡。
四）终止条件
有很多种终止条件的选择，各种不同的条件对算法的性能和解的质量有很大影响，本文只介绍一个常用的终止条件。即上一个最优解与最新的一个最优解的之差小于某个容差，即可停止此次马尔可夫链的迭代。
使用模拟退火法求函数 f(x,y) = 5sin(xy) + x2 + y2 的最小值
解：根据题意，我们设计冷却表进度表为：
即初始温度为 100
衰减参数为 0.95
马可夫链长度为 10000
Metropolis 的步长为 0.02
结束条件为根据上一个最优解与最新的一个最优解的之差小于某个容差。
使用METROPOLIS接受准则进行模拟, 程序如下
/*
* 模拟退火法求函数f(x,y) = 5sin(xy) + x^2 + y^2的最小值
* 结束条件为两次最优解之差小于某小量
*/
using System;
namespace SimulateAnnealing
{
class Class1
{
// 要求最优值的目标函数
static double ObjectFunction( double x, double y )
{
double z = 0.0;
z = 5.0 * Math.Sin(x*y) + x*x + y*y;
return z;
}
[STAThread]
static void Main(string[] args)
{
// 搜索的最大区间
const double XMAX = 4;
const double YMAX = 4;
// 冷却表参数
int MarkovLength = 10000; // 马可夫链长度
double DecayScale = 0.95; // 衰减参数
double StepFactor = 0.02; // 步长因子
double Temperature = 100; // 初始温度
double Tolerance = 1e-8; // 容差
double PreX,NextX; // prior and next value of x
double PreY,NextY; // prior and next value of y
double PreBestX, PreBestY; // 上一个最优解
double BestX,BestY; // 最终解
double AcceptPoints = 0.0; // Metropolis过程中总接受点
Random rnd = new Random();
// 随机选点
PreX = -XMAX * rnd.NextDouble() ;
PreY = -YMAX * rnd.NextDouble();
PreBestX = BestX = PreX;
PreBestY = BestY = PreY;
// 每迭代一次退火一次(降温), 直到满足迭代条件为止
do
{
Temperature *=DecayScale;
AcceptPoints = 0.0;
// 在当前温度T下迭代loop(即MARKOV链长度)次
for (int i=0;i<MarkovLength;i++)
{
// 1) 在此点附近随机选下一点
do
{
NextX = PreX + StepFactor*XMAX*(rnd.NextDouble()-0.5);
NextY = PreY + StepFactor*YMAX*(rnd.NextDouble()-0.5);
}
while ( !(NextX >= -XMAX && NextX <= XMAX && NextY >= -YMAX && NextY <= YMAX) );
// 2) 是否全局最优解
if (ObjectFunction(BestX,BestY) > ObjectFunction(NextX,NextY))
{
// 保留上一个最优解
PreBestX =BestX;
PreBestY = BestY;
// 此为新的最优解
BestX=NextX;
BestY=NextY;
}
// 3) Metropolis过程
if( ObjectFunction(PreX,PreY) - ObjectFunction(NextX,NextY) > 0 )
{
// 接受, 此处lastPoint即下一个迭代的点以新接受的点开始
PreX=NextX;
PreY=NextY;
AcceptPoints++;
}
else
{
double change = -1 * ( ObjectFunction(NextX,NextY) - ObjectFunction(PreX,PreY) ) / Temperature ;
if( Math.Exp(change) > rnd.NextDouble() )
{
PreX=NextX;
PreY=NextY;
AcceptPoints++;
}
// 不接受, 保存原解
}
}
Console.WriteLine("{0},{1},{2},{3}",PreX, PreY, ObjectFunction ( PreX, PreY ), Temperature);
} while( Math.Abs( ObjectFunction( BestX,BestY) – ObjectFunction (PreBestX, PreBestY)) > Tolerance );
Console.WriteLine("最小值在点:{0},{1}",BestX, BestY);
Console.WriteLine( "最小值为:{0}",ObjectFunction(BestX, BestY) );
}
}
}
l 结果：
最小值在点:-1.07678129318956,1.07669421564618
最小值为:-2.26401670947686
l 后记：
原来想写一系列的文章，介绍如何用C＃解数值问题，这是因为在CSDN上很少有数值计算方面的文章，所以希望能有所补充。
一开始在网上搜索模拟退火的资料并想作为 C ＃数值计算的一个例子，找不到现成的源码。后来自己实验了很久，终于将此程序写出，不敢私藏，拿出来作用模拟退火或者用 C ＃解数值算法问题的一个入门例子。
本文尽量避免太过学术化，如数学和物理名称和公式，仓促下笔，有很多地方可能讲得不是很清楚，希望各位体谅。任何问题或批评，可 EMAIL 与我： armylau2@163.com
另，模拟退火还可以应用到其它更多更复杂的问题，如“推销员问题”等组合优化问题。本例只是求一个二维函数的最小值问题，而且其冷却表参数的选择也过于简单，只能作用一个初步的入门简介，请读者注意。
327
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
