模拟退火算法是一种通用概率演算法，用来在一个大的搜索空间内找出最优解。相比于二分、三分等算法，模拟退火更加注意整体的情况，而非死磕在局部的变化。
爬山
可以拿爬山做例子：我们要找到山脉的最高峰，但是我们只知道眼前的点，哪边是下降的，但看不到远处的点是否上升。所以每次移动，我们随机选择一个方向。如果这个方向是上升的的（更优），那么就决定往那个方向走；如果这个方向是下降的（更差），那么“随机地接受”这个方向，接受就走，不接受就再随机一次。
模拟退火算法（流程）
- 随机产生一个初始解x0，令xbest＝ x0 ，并计算目标函数值E(x0);
- 设置初始温度T(0)=To，迭代次数i = 1;
- Do while T(i) > Tmin
① for j = 1~k
② 对当前最优解xbest按照某一邻域函数，产生一新的解xnew。计算新的目标函数值E(xnew) ，并计算目标函数值的增量ΔE = E(xnew) - E(xbest) 。
③ 如果ΔE ＜0，则xbest = xnew；
④ 如果ΔE ＞0，则p = exp(- ΔE /T(i));
（如果c = random[0,1] < p， xbest = xnew; 否则xbest = xbest）
⑤End for - i = i ＋ 1；
- End Do
-
输出当前最优点，计算结束。
Simulated Annealing
要注意的是，实际题目中exp(- ΔE /T(i))这个概率并不是必要情况，所以有时可以忽略（有点像贪心）
步骤：
①对于循环的判定可以设定为精度的判定，当步长STEP>EPS（EPS=1e-6）的时候执行循环，否则退出。
②在周围搜索出新的一个点（注意判定这个点是否满足规定范围）
③判定这个点是否为最优点
如果是则更新当前点
如果不是就以一定概率更新当前点（可忽略）
④缩小步长STEP
一、对于1维坐标寻找Y的最小值，可以先规定一个起点（这里是原点），和每次前进/后退的长度（STEP）。每次按照这个STEP前进/后退，一旦发现更优的点就更新。并且每次都按一定的比例缩小前进/后退的距离（降低再次搜索的范围）。当步数长度STEP小于规定精度时就可以退出搜索了。
const double EPS=1e-6;
const double r = 0.99;
const int dx[]= {-1,1};
…………………
double step=1;//规定初始步长长度
double x=0;//规定起始点
while(step > EPS)//设置精度范围
{
for(int i=0; i<2; i++)//分别前进和后退
{
double next_x = x+dx[i]*step;//确定下一个点的坐标x
if(F(next_x)<F(x))//判断是否为更优解
{
x=next_x;//更新坐标
}
}
step*=r;//降低再次搜索的范围
}
例题：http://acm.hdu.edu.cn/showproblem.php?pid=2899
F(x) = 6 * x7+8*x6+7x^3+5x^2-y*x (0 <= x <=100)
Can you find the minimum value when x is between 0 and 100.
输入Y值，求F（x）的最小值
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<ctime>
using namespace std;
const double EPS=1e-6;
const double r = 0.99;
const int dx[]= {-1,1};
double y;
double F(double x)
{
return 6*pow(x,7)+8*pow(x,6)+7*pow(x,3)+5*pow(x,2)-y*x;
}
int main()
{
int T;
scanf("%d", &T);
while(T--)
{
scanf("%lf", &y);
double step=1;
double x=0;
while(step > EPS)
{
for(int i=0; i<2; i++)
{
double next_x = x+dx[i]*step;
if(F(next_x)<F(x))
{
x=next_x;
}
}
step*=r;
}
printf("%.4lf\n", F(x));
}
return 0;
}
二、对于2维坐标寻找Z的最小值
同样可以先规定一个起点（这里还是原点），和每次行进的长度（STEP）。每次按照这个STEP向东、南、西、北、东南、西南、东北、西北八个方向行进，一旦发现更优的点就更新。并且每次都按一定的比例缩小行进的距离（降低再次搜索的范围）。当步数长度STEP小于规定精度时就可以退出搜索了。
const double EPS=1e-6;
const int dx[]= {-1,-1,-1,0,0,1,1,1};
const int dy[]= {-1,0,1,-1,1,-1,0,1};
const double r = 0.99;
……………………
double step=1.0;//规定初始步长长度
double x=0, y=0;//规定初始点
double z=F(x, y);
while(step > EPS)
{
for(int i=0; i<8; i++)//向8个方向分别遍历
{
double next_x = x + dx[i] * step;
double next_y = y + dy[i] * step;
double next_z = F(next_x, next_y)
//这里还要判断next_z是否满足规定区域
if(distance(next_x, next_y, next_z) < distance(x, y, z))//判断更新最优解的条件
{
x = next_x;//满足条件就更新坐标
y = next_y;//满足条件就更新坐标
z = next_z;//满足条件就更新坐标
}
}
step*=r;//降低再次搜索的范围
}
例题:http://acm.hdu.edu.cn/showproblem.php?pid=5017
给定椭球公式，求从原点到椭球的最短距离
椭球公式
#include<cstdio>
#include<cmath>
using namespace std;
const double EPS=1e-6;
const double INF=0x3f3f3f3f;
const int dx[]= {-1,-1,-1,0,0,1,1,1};
const int dy[]= {-1,0,1,-1,1,-1,0,1};
const double r = 0.99;
double a, b, c, d, e, f, x, y, z;
double dis(double x, double y, double z)
{
return sqrt(x*x+y*y+z*z);
}
double F(double x, double y)
{
double aa, bb, cc;
aa=c;
bb=d*y+e*x;
cc=a*x*x+b*y*y+f*x*y-1;
if(bb*bb-4*aa*cc<0) return -INF;//点Z不在椭球面上
if( dis(x,y,(sqrt(bb*bb-4*aa*cc)-bb))/(2*aa)>dis(x,y,(-sqrt(bb*bb-4*aa*cc)-bb))/(2*aa))//判断一元二次方程到底取哪个解
return (-sqrt(bb*bb-4*aa*cc)-bb)/(2*aa);
else
return (sqrt(bb*bb-4*aa*cc)-bb)/(2*aa);
}
int main()
{
while(~scanf("%lf%lf%lf%lf%lf%lf", &a, &b, &c, &d, &e, &f))
{
double step=1.0;
double x=0, y=0;
double z=F(x, y);
while(step > EPS)
{
for(int i=0; i<8; i++)
{
double next_x = x + dx[i] * step;
double next_y = y + dy[i] * step;
double next_z = F(next_x, next_y);
if(next_z >= INF) continue;//点Z不在椭球面上
if(dis(next_x, next_y, next_z) < dis(x, y, z))
{
x = next_x;
y = next_y;
z = next_z;
}
}
step*=r;
}
printf("%lf\n", dis(x, y, z));
}
return 0;
}
当然在一定情况下，上面的搜索并不能满足要求，首先是对于初始点的位置，一个初始点可能会限定搜索的范围，最后无法找出最优解。其次是方向的选择，8个方向也并不是最优的方向。对于改进方式，就是在规定的区域内，随机分布多个点，并且在每个点附近寻找最优解的时候，采取任意方向搜索。
三、对于二维规定区域找出满足条件的最优解（坐标）
思路和上面大致相同，只是将原来的单点逐步搜索改为多点逐步搜索，搜索方向也变为任意方向搜索。
规定区域内找随机分布的点
for(i=0; i<num; i++)
{
m[i].x=(rand()%1000+1)/1000.0*x0;
m[i].y=(rand()%1000+1)/1000.0*y0;
}
初始步长的长度
double step=sqrt(x0*x0+y0*y0)/2;
任意方向的处理
double angle=(rand()%1000+1)/1000.0*2*PI;
还要注意判定每次新发现的点（最优解）是否在规定区域内部喔~
例题：http://acm.hdu.edu.cn/showproblem.php?pid=3932
题目大意：在平面中的已知点中，找到一个点A，这个点要求是到其他所有最长点的最短情况。（求出距离最远的那个点到A的长度,并且这个长度是情况中最短的长度）输出A的坐标和A到最远点的距离。
#include<cstdio>
#include<ctime>
#include<cmath>
#include<algorithm>
using namespace std;
const double EPS=1e-6;
const double PI=acos(-1);
const double INF=0x3f3f3f3f;
const double r=0.8;
typedef struct st
{
double x, y;
} ST;
ST a[10005];
ST m[10005];
int n;
double d[10005];
double dis(ST A, ST B)
{
return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y));
}
double MAX(ST t)
{
int i;
double maxx=0;
for(i=0; i<n; i++)
{
maxx=max(maxx,dis(t,a[i]));
}
return maxx;
}
int main()
{
double x0, y0;
srand(time(NULL));
int init_num=15;
while(~scanf("%lf%lf%d", &x0, &y0, &n))
{
int i, j;
for(i=0; i<n ;i++)
{
scanf("%lf%lf", &a[i].x, &a[i].y);
}
for(i=0; i<init_num; i++)
{
m[i].x=(rand()%1000+1)/1000.0*x0;
m[i].y=(rand()%1000+1)/1000.0*y0;
d[i]=MAX(m[i]);
}
double step=sqrt(x0*x0+y0*y0)/2;
ST next, temp;
while(step>EPS)
{
for(i=0; i<init_num; i++)
{
temp.x=m[i].x;
temp.y=m[i].y;
for(j=0; j<50; j++)
{
double angle=(rand()%1000+1)/1000.0*2*PI;
next.x=temp.x+cos(angle)*step;
next.y=temp.y+sin(angle)*step;
if(next.x<0 || next.x>x0 || next.y<0 || next.y>y0) continue;
if(MAX(next)<d[i])
{
m[i].x=next.x;
m[i].y=next.y;
d[i]=MAX(next);
//printf("ok\n");
}
}
}
step*=r;
}
int flag;
double ans=INF;
for(i=0; i<init_num; i++)
if(d[i]<ans)
{
ans=d[i];
flag=i;
}
printf("(%.1lf,%.1lf).\n",m[flag].x,m[flag].y);
printf("%.1lf\n",ans);
}
return 0;
}
网友评论