C#数值计算之模拟退火法简介(二)

在上一篇文章中讲述了模拟退火的基本原理，以下以一个实际的例子来说明，其中所有的源码已贴出，可以从中了解到很多细节。

使用模拟退火法求函数f(x,y) = 5sin(xy) + x2 + y2的最小值

解：根据题意，我们设计冷却表进度表为：

即初始温度为100

衰减参数为0.95

马可夫链长度为10000

Metropolis的步长为0.02

结束条件为根据上一个最优解与最新的一个最优解的之差小于某个容差。

使用METROPOLIS接受准则进行模拟, 程序如下

/*

* 模拟退火法求函数f(x,y) = 5sin(xy) + x^2 + y^2的最小值

* 日期:2004-4-16

* 作者:ARMYLAU

*　EMAIL:armylau2@163.com

* 结束条件为两次最优解之差小于某小量

*/

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

另，模拟退火还可以应用到其它更多更复杂的问题，如“推销员问题”等组合优化问题。本例只是求一个二维函数的最小值问题，而且其冷却表参数的选择也过于简单，只能作用一个初步的入门简介，请读者注意。

l 参考文献：

1． http://www.computer-dictionary-online.org/index.asp?q=simulated+annealing　计算机词典
2． Numeric Recipes in C

3． 计算方法丛书 非数值并行算法 （第一册） 模拟退火算法