[HDU1109]模拟退火算法

模拟退火的基本思想:

(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)，每个T值的迭代次数L

(2) 对k=1，……，L做第(3)至第6步：

(3) 产生新解$Sprime $

(4) 计算增量$Delta tprime  = Cleft( {Sprime } ight) - Cleft( S ight)$，其中$Cleft( S ight)$为评价函数

(5) 若$Delta tprime  < 0$则接受$Sprime $作为新的当前解，否则以概率$expleft( { - Delta tprime /T} ight)$接受$Sprime $作为新的当前解.

(6) 如果满足终止条件则输出当前解作为最优解，结束程序。 终止条件通常取为连续若干个新解都没有被接受时终止算法。

(7) T逐渐减少，且T->0，然后转第2步。

从算法的流程上看，模拟退火算法包括三函数、两准则，即状态产生函数、状态接受函数、温度更新函数、内循环终止准则和外循环终止准则，这些环节的设计将决定模拟退火算法的优化性能。此外，初温的选择对模拟退火算法性能也有很大的影响。

状态产生函数：

   原则：设计状态产生函数(领域函数)的出发点应该是尽可能保证产生的候选解遍布全部的解空间。通常，状态产生函数由两部分组成，即产生候选解的方式和候选解吃剩的概率分布。
   方法：在当前状态下的领域结构内以一定的概率方式(均匀分布、正态分布、指数分布等)产生

状态接受函数：

   原则：函数一般以概率的方式给出，不同接受函数的差别主要在于接受概率的形式不同。设计状态接受概率，应该遵循以下原则：

1）在固定温度下，接受使目标函数下降的候选解的概率要大于使目标函数上升的候选解的概率； 
2）随温度的下降，接受目标函数上升的解的概率要逐渐减少； 
3）当温度趋于零时，只能接受目标函数下降的解。方法：状态接受函数的引入是模拟退火算法实现全局搜索的最关键的因素，模拟退火算法中通常用作为状态接受函数。

初始温度、温度更新函数、内循环终止准则和外循环终止准则通常被称为退火历程。

初始温度：

    原则：通过理论分析可以得到初温的解析式，但解决实际问题难以得到精确的参数；实际应用往往要让初温充分大。实验表明：初温越大，获得高质量解的机率越大，但花费较多的计算时间。

方法： 
（1）均匀抽样一组状态，以各状态目标值的方差为初温； 
（2）随机产生一组状态，确定两两状态间的最大目标差值，根据差值，利用一定的函数确定初温，譬如 ,其中为初始接受函数； 
（3）利用经验公式。

温度更新函数：

 温度更新函数，即温度下降方式，用于在啊外循环中修改温度值。

内循环终止准则：

常用的Metropolis抽样准则： 
1）检验目标函数的均值是否稳定； 
2）连续若干步的目标值变化较小； 
3）按一定步数抽样。

外循环终止准则：

 1）设置终止温度的阈值；
 2）设置外循环迭代次数；
 3）算法搜索到的最优值连续若干步保持不变；
 4）概率分析方法。

[p(dE) = exp(dE/kT)]其中k是一个常数，exp表示自然指数，且dE<0。这条公式说白了就是：温度越高，出现一次能量差为dE的降温的概率就越大；温度越低，则出现降温的概率就越小。又由于dE总是小于0（否则就不叫退火了），因此dE/kT < 0 ，所以P(dE)的函数取值范围是(0,1) 。

 若ΔT<0则接受S′作为新的当前解S，否则以概率exp(-ΔT/T)接受S′作为新的当前解S。

http://acm.hdu.edu.cn/showproblem.php?pid=1109

题意：给定X,Y,M,求区域内点到已知点最短距离的最大值。

解题关键：

1、注意一定不要忘记设立初始解。

2、T_min、T_max转化为步长，以保证精确度。

3、内循环、外循环一定要理解。

4、此程序其实还有一个问题未解决，即以一定的概率接受非最优解，而由于未找到良好的评估函数，故可以直接用时间之比即可，若有时候WA时，尝试一下下面第二个程序。

5、在这里其实就是将外循环的降温当做误差范围。

6、第一个程序相当于多点爬山？不算随机算法。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<cmath>
#include<iostream>
using namespace std;
typedef long long ll;
#define INF 100000000
#define N 30
#define PI acos(-1.0)
double X,Y,M;
double step,eps=1e-3,k=0.55;//step代表最长的跨度 ,注意这种用法 ，用来控制精度
struct point{
    double x,y;
}people[10002],num[10002];
double d[10002];

double dist(point a,point b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

bool judge(point a){
    return a.x>=0&&a.x<=X&&a.y>=0&&a.y<=Y;
}

double rand1(double l,double r){//随机数生成函数
    return rand()%10000/10000.0*(r-l)+l;
}

double fun(point a){
    double td=INF;
    for(int i=0;i<M;i++){
        td=min(td,dist(num[i],a));
    }
    return td;
}
int main(){
    int t;
    scanf("%d",&t);
    while(t--){
        double max1=-1.0;
        int maxn=0;
        cin>>X>>Y>>M;
        for(int i=0;i<M;i++){
            scanf("%lf%lf",&num[i].x,&num[i].y);
        }
        
        for(int i=0;i<N;i++){
            people[i].x=rand1(0,(double)X);
            people[i].y=rand1(0,(double)Y);
            d[i]=fun(people[i]);
        }
        
        step=max(X,Y);
        while(step>eps){
            for(int i=0;i<N;i++){//初始状态一般为30个
                for(int j=0;j<50;j++){//一般循环50次左右，看题目要求的时间可以变化
                    point tem;
                    double angle=rand1(0,2*PI);//枚举任何角度，使得到的新点向四周扩散
                    tem.x=people[i].x+cos(angle)*step;
                    tem.y=people[i].y+sin(angle)*step;
                    if(!judge(tem)) continue;
                    double dd=fun(tem);
                    if(dd>d[i]){     //这里一定要注意，必须要更新，不要忘记
                        d[i]=dd;
                        people[i]=tem;
                    }
                }
            }
            step*=k;
        }
        for(int i=0;i<N;i++){//找到退火后的状态中，最优的解
            if((d[i]-max1)>eps){
                maxn=i;
                max1=d[i];
            }
        }
        printf("The safest point is (%.1f, %.1f).
",people[maxn].x,people[maxn].y);
    }
}

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define INF 100000000
#define EPS 1e-6
#define N 30
#define PI acos(-1.0)
double X,Y,M;
double step,eps=1e-3,k=0.55;//step代表最长的跨度 ,注意这种用法 ，用来控制精度 
struct point{
    double x,y;
}people[10002],num[10002];
double d[10002];

double dist(point a,point b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

bool judge(point a){
    return a.x>=0&&a.x<=X&&a.y>=0&&a.y<=Y;
}

double rand1(double l,double r){//随机数生成函数 
    return rand()%10000/10000.0*(r-l)+l;
}

double fun(point a){
    double td=INF;
    for(int i=0;i<M;i++){
        td=min(td,dist(num[i],a));
    }
    return td;
}
int main(){
    int t;
    scanf("%d",&t);
    while(t--){
        double max1=-1.0;
        int maxn=0;
        cin>>X>>Y>>M;
        for(int i=0;i<M;i++){
            scanf("%lf%lf",&num[i].x,&num[i].y);
        }
        
        for(int i=0;i<N;i++){
            people[i].x=rand1(0,(double)X);
            people[i].y=rand1(0,(double)Y);
            d[i]=fun(people[i]);
        }
        
        step=max(X,Y);
        while(step>eps){
            for(int i=0;i<N;i++){//初始状态一般为30个 
                for(int j=0;j<50;j++){//一般循环50次左右，看题目要求的时间可以变化 
                    point tem;
                    double angle=rand1(0,2*PI);//枚举任何角度，使得到的新点向四周扩散 
                    tem.x=people[i].x+cos(angle)*step;
                    tem.y=people[i].y+sin(angle)*step;
                    if(!judge(tem)) continue;
                    double dd=fun(tem);
                    if(dd>d[i]){     //这里一定要注意，必须要更新，不要忘记 
                        d[i]=dd;
                        people[i]=tem;
                    }
                    else{
                        if(rand()%10000/10000.0>exp(k-1)){
                            d[i]=dd;
                            people[i]=tem;
                        }
                    }
                } 
            }
            step*=k;
        }
        for(int i=0;i<N;i++){//找到退火后的状态中，最优的解 
            if((d[i]-max1)>EPS){
                maxn=i;
                max1=d[i];
            }
        }
        printf("The safest point is (%.1f, %.1f).
",people[maxn].x,people[maxn].y);
    }
}