模拟退火算法求函数最小值问题的C/C++实现

1、模拟退火算法简介

模拟退火算法(Simulated Annealing，SA)最早的思想是由N. Metropolis等人于1953年提出。1983 年,S. Kirkpatrick 等成功地将退火思想引入到组合优化领域。它是基于Monte-Carlo迭代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是一种通用的优化算法，理论上算法具有概率的全局优化性能,目前已在工程中得到了广泛应用，诸如VLSI、生产调度、控制工程、机器学习、神经网络、信号处理等领域。

2、模拟退火算法的基本原理

模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis准则，粒子在温度T时趋于平衡的概率为e(-ΔE/(kT))，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件：最终温度。

3、模拟退火算法的伪代码

(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)， 每个T值的迭代次数L 
(2) 对k=1，……，L做第(3)至第6步： 
  (3) 产生新解S′ 
  (4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数 
  (5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/T)接受S′作为新的当前解. 
  (6) 如果满足终止条件则输出当前解作为最优解，结束程序。 
        终止条件通常取为连续若干个新解都没有被接受时终止算法。 
(7) T逐渐减少，且T->最终温度，然后转第2步。

4、应用模拟退火算法求解函数最小值C/C++实现如下：

#include <iostream>
#include <math.h>
#include <time.h>
//f(x)=11*sin(6*x)+7*cos(5*x),x∈[0,2*pi]，求最小值
#define pi 3.14159
#define num 30000 //迭代次数
double k=0.1;
double r=0.9; //用于控制降温的快慢
double T=2000; //系统的温度，系统初始应该要处于一个高温的状态
double T_min =20;//温度的下限，若温度T达到T_min，则停止搜索
//返回指定范围内的随机浮点数
double rnd(double dbLow,double dbUpper)
{
    double dbTemp=rand()/((double)RAND_MAX+1.0);
    return dbLow+dbTemp*(dbUpper-dbLow);
}
double func(double x)//目标函数
{
	return 7*sin(6*x)+6*cos(5*x);
}
int main()
{	
	double best=func(rnd(0.0,2*pi));
	double dE,current;	
	int i;
	while( T > T_min )
    {
		for(i=0;i<num;i++)
		{
		//用当前时间点初始化随机种子，防止每次运行的结果都相同
		time_t tm;
		time(&tm);
		unsigned int nSeed=(unsigned int)tm;
		srand(nSeed);

		current=func(rnd(0.0,2*pi));
		dE = current - best ;

		if ( dE < 0 ) //表达移动后得到更优解，则总是接受移动
			best = current ; 
		else
		{
			// 函数exp( dE/T )的取值范围是(0,1) ，dE/T越大，则exp( dE/T )也越大
			if ( exp( dE/(T*k) ) > rnd( 0.0 , 1.0 ) )
				best = current ; 
		}
		}		
		T = r * T ;//降温退火 ，0<r<1 。r越大，降温越慢；r越小，降温越快
    }
	printf("最小值是 %f\n",best);
	return 0;
}

运行结果：