题解 P1337 【[JSOI2004]平衡点 / 吊打XXX】

· · 题解

关于模拟退火

简介

模拟退火是一种随机化算法。
对于一个当前最优解附近的非最优解，爬山算法直接舍去了这个解。而很多情况下，我们需要去接受这个非最优解从而跳出这个局部最优解，即为模拟退火算法。
当一个问题的方案数量极大（甚至是无穷的）而且不是一个单峰函数时，常使用模拟退火求解。

实现

如果新状态的解更优则修改答案，否则以一定概率接受新状态。

模拟退火时有三个参数：初始温度$T_0$，降温系数$d$ ，终止温度$T_k$。 $T_0$是一个比较大的数（至少需要大于搜索范围）
$d$ 是一个非常接近1但是小于1的数
$T_k$是一个接近0的正数
定义当前温度为$T$，新状态与已知状态（由已知状态通过随机的方式得到）之间的能量差为$\Delta E$，（当要求的$E$为最大值时）发生状态转移（修改最优解）的概率为
$P(\Delta E)=\begin{cases}1,\quad \Delta E>0\\e^{\frac{\Delta E}{T}},\quad \Delta E\leq 0\end{cases}$
首先让温度$T=T_0$，进行转移尝试，再让$T=T*d$。当$T<T_k$时模拟退火过程结束，当前最优解即为最终的最优解。
有时为了使得到的解更有质量，会在模拟退火结束后，以当前温度在得到的解附近多次随机状态，尝试得到更优的解（其过程与模拟退火相似）。

如图随着温度的降低，跳跃越来越不随机，最优解也越来越稳定

关于本题

虽然模拟退火不是正解，但是可以AC!
而且本题也是一道非常优秀的模拟退火例题!

当重力势能最小时稳定，将物理模型简化
设桌子的高度为$h$，绳子的长度为$l$,每个点到$x$的距离为$dis[i]$
重力势能为$\sum_{i=1}^{n}(h-(l-dis[i]))*W[i]=\sum_{i=1}^{n}(h-l+dis[i])*W[i]$
由于$h-l$为定值，所以重力势能最小的情况等价于$\sum_{i=1}^{n}dis[i]*W[i]$最小
即 求$n$个点的带权费马点

Code:
#include<bits/stdc++.h>
#define For(i,a,b) for(int i=(a),i##_=(b);i<=i##_;++i)
#define M 1005
#define db double
using namespace std;
int X[M],Y[M],W[M],n;
db ans,ansx,ansy;
db Pow(db x){
    return x*x;
}
db Dis(db x,db y){
    return sqrt(Pow(x)+Pow(y));
}
db Rand(){
    return (db)rand()/RAND_MAX;//返回一个[0,1]之前的小数
}
db Cal(db x,db y){
    db res=0;
    For(i,1,n)res+=Dis(x-X[i],y-Y[i])*W[i];
    //在Cal中更新答案
    if(res<ans)ans=res,ansx=x,ansy=y;
    return res;
}
void SimulateAnneal(){
    db T0=100000,Tk=0.001,d=0.97;
    db nowx=ansx,nowy=ansy,T=T0;
    while(T>Tk){
        //Rand()*2-1 返回一个[-1,1]之间的数
        db nxtx=nowx+T*(Rand()*2-1);
        db nxty=nowy+T*(Rand()*2-1);
        db tmp=Cal(nxtx,nxty)-Cal(nowx,nowy);
        //tmp<0 即新解更优，那么一定接受新解
        //否则 以一定概率接受新解
        if(tmp<0||exp(-tmp/T)>Rand())nowx=nxtx,nowy=nxty;
        T*=d;
    }
    //以当前温度在得到的解附近多次随机状态，尝试得到更优的解
    int t=1000;
    For(i,1,t){
        db nxtx=ansx+T*(Rand()*2-1);
        db nxty=ansy+T*(Rand()*2-1);
        Cal(nxtx,nxty);
    }
}
int main(){
    srand(19260817);
    scanf("%d",&n);
    For(i,1,n){
        scanf("%d%d%d",&X[i],&Y[i],&W[i]);
        ansx+=X[i],ansy+=Y[i];
    }
    //初始值可以随机，在本题中采用所有点坐标的平均值
    ansx/=n,ansy/=n;
    ans=Cal(ansx,ansy);
    SimulateAnneal();
    printf("%.3f %.3f\n",ansx,ansy);
    return 0;
}
//这份代码不保证一次AC哦！
//可以手调参数来加深理解哦！