模拟退火算法介绍(摘自08集训队 顾研《浅谈随机化思想在几何问题中的应用》)


    一.模拟退火算法的原理

模拟退火算法是一种元启发式(Meta-Heuristics)算法,来源于固体退火原理,将固体加热至充分高的温度,再让其徐徐冷却。加热时,固体内部粒子随温升变为无序状,内能增大,而徐徐冷却时粒子渐趋有序,在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小。根据Metropolis准则,粒子在温度T时趋于平衡的概率为e-ΔE/kt,其中E为温度T时的内能,ΔE为其改变量,kBoltzmann常数。

    二.模拟退火算法的模型

① 初始化:初始温度T(充分大),初始解状态S(算法迭代的起点), 每次迭代次数L

② for k=1 to L 做③至⑥

③ 产生新解S’

④ 计算增量Δt′=C(S′)-C(S),其中C(S)为评价函数

⑤ 若Δt′<0则接受S’作为新的当前解,否则以概率e-Δt/k接受S’作为新的当前解

⑥ 如果满足终止条件则输出当前解作为最优解,结束程序

⑦ T逐渐减少,然后转②

    回到此题,题意为:平面上有一个矩形,在矩形内有一些陷阱。求得矩形内一个点,该点离与它最近的已知陷阱最远(点的个数1000)。精度要求:1/10

这题的精度要求不高,一个很朴素的想法是枚举平面中的一些网格点并进行判断。当然这样的点太多了,我们必须选一些比较有“希望”的点,同时还不能忽略任何平面上的任何一点。

我们使用类比的方法引入模拟退火算法,初始解状态S可以在矩形内随便选取,初始温度T对应于每次调整的距离D,产生新解的方式是在目前解为圆心、半径为D的圆周上任取一点,评价函数取距离最近的陷阱的距离,终止条件为D足够小。

由此可得本问题的模拟退火算法:由初始解S和控制参数初值D开始,对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代,并逐步衰减D值,算法终止时可以得到一组解,这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制,包括控制参数的初值D及其衰减因子D、每个值D时的迭代次数L

模拟退火算法还有一个特点:具有并行性。因此我们可以将初始解S改成初始解集,对于每个候选解都进行迭代,答案取最终解集的最优解。

由于我们必须保证能覆盖矩形上的每个位置,因此在①中确定p后(我们按网格状放置),②中的delta可取max{矩形边长}/sqrt(n)。设算法的迭代次数为t,则算法的复杂度O(P*L*t*n)
    


#include<iostream>
#include
<time.h>
#include
<cmath>
#define sqr(x) (x)*(x)
#define FOR(i,a,b) for(int i=a;i<b;i++)
#define FF(i,a) for(int i=0;i<a;i++)
using namespace std;
const double eps=1e-3;
const double INF=1e20;
const double pi=acos(-1.0);
struct Houxuan
{
    
double x,y,dist;
}pp[
30];
struct PT
{
    
double x,y;
}hole[
1000];
double dis(double x1,double y1,double x2,double y2)
{
    
return sqrt(sqr(x1-x2)+sqr(y1-y2));
}
int n,cas;
int X,Y;
int main()
{
    
int P=15,L=25;
    scanf(
"%d",&cas);
    srand((unsigned)time(NULL));
    
while(cas--)
    {
        scanf(
"%d%d%d",&X,&Y,&n);
        
double delta=double(max(X,Y))/(sqrt(1.0*n));
        FF(i,n)
            scanf(
"%lf%lf",&hole[i].x,&hole[i].y);
        FF(i,P)
        {
            pp[i].x
=double(rand()%1000+1)/1000.000*X;
            pp[i].y
=double(rand()%1000+1)/1000.000*Y;
            pp[i].dist
=INF;
            FF(j,n)
                pp[i].dist
=min(pp[i].dist,dis(pp[i].x,pp[i].y,hole[j].x,hole[j].y));
        }
        
while(delta>eps)
        {
            FF(i,P)
            {
                Houxuan tmp
=pp[i];
                FF(j,L)
                {
                    
double theta=double(rand()%1000+1)/1000.000*10*pi;
                    Houxuan t;
                    t.x
=tmp.x+delta*cos(theta);
                    t.y
=tmp.y+delta*sin(theta);
                    
if(t.x>X||t.x<0||t.y>Y||t.y<0)
                        
continue;
                    t.dist
=INF;
                    FF(k,n)
                    {
                        t.dist
=min(t.dist,dis(t.x,t.y,hole[k].x,hole[k].y));
                    }
                    
if(t.dist>pp[i].dist)
                        pp[i]
=t;                    
                }
            }
            delta
*=0.8;
        }
        
int idx=0;
        FOR(i,
1,P)
        {
            
if(pp[i].dist>pp[idx].dist)
                idx
=i;
        }    
        printf(
"The safest point is (%.1lf, %.1lf).\n",pp[idx].x,pp[idx].y);    
    }
    
return 0;
}