膜模魔 模拟退火

一. 爬山算法 ( Hill Climbing )

介绍模拟退火前，先介绍爬山算法。爬山算法是一种简单的贪心搜索算法，该算法每次从当前解的临近解空间中选择一个最优解作为当前解，直到达到一个局部最优解。

爬山算法实现很简单，其主要缺点是会陷入局部最优解，而不一定能搜索到全局最优解。如图1所示：假设C点为当前解，爬山算法搜索到A点这个局部最优解就会停止搜索，因为在A点无论向那个方向小幅度移动都不能得到更优的解。

二. 模拟退火(SA,Simulated Annealing)思想

爬山法是完完全全的贪心法，每次都鼠目寸光的选择一个当前最优解，因此只能搜索到局部的最优值。模拟退火其实也是一种贪心算法，但是它的搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。以图1为例，模拟退火算法在搜索到局部最优解A后，会以一定的概率接受到E的移动。也许经过几次这样的不是局部最优的移动后会到达D点，于是就跳出了局部最大值A。

模拟退火算法描述：

若J( Y(i+1) )>= J( Y(i) )  (即移动后得到更优解)，则总是接受该移动

若J( Y(i+1) )< J( Y(i) )  (即移动后的解比当前解要差)，则以一定的概率接受移动，而且这个概率随着时间推移逐渐降低（逐渐降低才能趋向稳定）

这里的“一定的概率”的计算参考了金属冶炼的退火过程，这也是模拟退火算法名称的由来。

根据热力学的原理，在温度为T时，出现能量差为dE的降温的概率为P(dE)，表示为：

P(dE) = exp( dE/(kT) )

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

随着温度T的降低，P(dE)会逐渐降低。

我们将一次向较差解的移动看做一次温度跳变过程，我们以概率P(dE)来接受这样的移动。

下面给出模拟退火的伪代码表示。

/*
* J(y)：在状态y时的评价函数值
* Y(i)：表示当前状态
* Y(i+1)：表示新的状态
* r： 用于控制降温的快慢
* T： 系统的温度，系统初始应该要处于一个高温的状态
* T_min ：温度的下限，若温度T达到T_min，则停止搜索
*/
while( T > T_min )
{
　　dE = J( Y(i+1) ) - J( Y(i) ) ; 
 
　　if ( dE >=0 ) //表达移动后得到更优解，则总是接受移动
       　　Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
　　else{    // 函数exp( dE/T )的取值范围是(0,1) ，dE/T越大，则exp( dE/T )也
　　　　if ( exp( dE/T ) > random( 0 , 1 ) )
　　　　Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
　　}
　　T = r * T ; //降温退火 ，0<r<1 。r越大，降温越慢；r越小，降温越快
　　/*
　　* 若r过大，则搜索到全局最优解的可能会较高，但搜索的过程也就较长。若r过小，则搜索的过程会很快，但最终可能会达到一个局部最优值
　　*/
　　i ++ ;
}

注意事项：

1. 温度T的初始值设置问题。 初始温度高，则搜索到全局最优解的可能性大，但因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。
2. 退火速度r问题。 模拟退火算法的全局搜索性能也与退火速度密切相关。同一温度下的“充分”搜索(退火)是相当必要的，但这需要计算时间。
3. 温度管理问题 降温系数应为正的略小于1.00的常数。

有几题裸题：

BZOJ 3680

这道题，一切自然变化进行的方向都是使能量降低，因为能量较低的状态比较稳定。

因为物重一定，绳子越短，重物越低，势能越小，势能与物重成正比，所以，只要使得也就是总的重力势能最小，就可以使系统平衡

#include<cstdio>
#include<cctype>
#include<cmath>
#include<cstdlib>

struct Node
{
    int x,y,weight;
}node[10005];
int n;

double potential_energy(double nowx,double nowy)
{
    double sum=0;
    for(int i=1;i<=n;i++)
    {
        double delx=nowx-node[i].x;
        double dely=nowy-node[i].y;
        sum+=(sqrt(delx*delx+dely*dely))*node[i].weight;
        //物重一定，绳子越短，重物越低，势能越小 
        //势能又与物重成正比 
    }
    return sum;//在(nowx,nowy)处的总势能 
}

double xans,yans;//最终答案 
double ans=1e18+7,t;//势能与温度
const double delta=0.993;//降温系数

void simulate_anneal()
{
    double xx=xans;//钦定一个初始位置
    double yy=yans;
    t=1926;//t是温度
    while(t>1e-14)
    {
        double xtemp=xans+(rand()*2-RAND_MAX)*t;
        double ytemp=yans+(rand()*2-RAND_MAX)*t;
        //随机一个新的坐标，变化幅度为t
        //这里要注意rand()和rand()*2-RAND_MAX的区别
        //rand()的范围是0~RAND_MAX-1
        //rand()*2-RAND_MAX的范围是-RAND_MAX到RAND_MAX-1 

        double new_ans=potential_energy(xtemp,ytemp);//计算当前解的势能 
        double DE=new_ans-ans;
        if(DE<0)//如果是一个更优解
        {
            xx=xtemp;
            yy=ytemp;//就接受
            xans=xx;
            yans=yy;
            ans=new_ans;
        }
        else if(exp(-DE/t)*RAND_MAX>rand())//能否接受这个差 
        {
            //更新坐标
            xx=xtemp;
            yy=ytemp;
        }
        t*=delta;//降温 
    }
}

void SA()//洗把脸就AC了 
{
    simulate_anneal();
    simulate_anneal();
    simulate_anneal();
}

int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        scanf("%lf%lf%lf",&node[i].x,&node[i].y,&node[i].weight);
        node[i].x=Read();
        node[i].y=Read();
        node[i].weight=Read();
    }
    SA();
    printf("%.3lf %.3lf",xans,yans);
    return 0;
}

POJ2680

这题算是简化版的模拟退火，并没有用到：以概率P(dE)来接受这样的移动。

所以一开始用普通的模拟退火来做卡了半天，最后发现其实是简单的移动。

我们首先任选一个点作为球心，并找到点集中与它距离最远的点，我们让球心靠近最远的点，不断重复此过程，就可以让球心达到稳定态了！此时我们就找到了最小球。

代码：

#include <iostream>  
#include<cstdio>  
#include<algorithm>  
#include<cstring>  
#include<cmath>  
using namespace std;  
const double eps=1e-7;  
struct point3D  
{  
    double x,y,z;  
} data[35];  
int n;  
double dis(point3D a,point3D b)  
{  
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z));  
}  
double solve()  
{  
    double step=100,ans=1e30,mt;  
    point3D z;  
    z.x=z.y=z.z=0;  
    int s=0;  
    while(step>eps)  
    {  
        for(int i=0; i<n; i++)  
            if(dis(z,data[s])<dis(z,data[i])) s=i;  
        mt=dis(z,data[s]);  
        ans=min(ans,mt);  
        z.x+=(data[s].x-z.x)/mt*step;  
        z.y+=(data[s].y-z.y)/mt*step;  
        z.z+=(data[s].z-z.z)/mt*step;  
        step*=0.98;  
    }  
    return ans;  
}  
int main()  
{ // freopen("t.txt","r",stdin);
    double ans;  
    while(~scanf("%d",&n),n)  
    {  
        for(int i=0; i<n; i++)  
            scanf("%lf%lf%lf",&data[i].x,&data[i].y,&data[i].z);  
        ans=solve();  
        printf("%.5f\n",ans);  
    }  
    return 0;  
}