C#數值計算之類比退火法簡介(二)
來源：互聯網
上載者：User
在上一篇文章中講述了類比退火的基本原理，以下以一個實際的例子來說明，其中所有的源碼已貼出，可以從中瞭解到很多細節。
使用類比退火法求函數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． 計算方法叢書 非數值並行演算法 （第一冊） 類比退火演算法
