在上一篇文章中講述了模擬退火的基本原理，以下以一個實際的例子來說明，其中所有的源碼已貼出，可以從中了解到很多細節。
使用模擬退火法求函數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:[email protected]
* 結束條件為兩次最優解之差小于某小量
*/
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與我：[email protected]
另，模擬退火還可以應用到其它更多更復雜的問題，如“推銷員問題”等組合優化問題。本例只是求一個二維函數的最小值問題，而且其冷卻表參數的選擇也過于簡單，只能作用一個初步的入門簡介，請讀者注意。
l 參考文獻：
1． http://www.computer-dictionary-online.org/index.asp?q=simulated+annealing 計算機詞典
2． numeric recipes in c
3． 計算方法叢書 非數值并行算法 （第一冊） 模擬退火算法
菜鳥學堂：
