基于模拟退火的随机性算法的应用&广义费马点丨BZOJ3680

模拟退火

总的来说只需要知道它是一种启发式搜索算法，并且能以很大概率获得全局最优解。

我们来简单说一下这个算法：在铁匠打铁时，需要将铁加热至高温，将铁水加至容器后，降温至室温，铁原子便会排列有序，形成一个铁具。降温的过程叫做退火过程，这个物理过程 中能量是从功能量变为低能量的。

而科学家们发现，若在当前温度T下，能量为E(0)，在降温单位时间后，能量$E(1)$并不一定比$E(0)$小，$E(1)=x$这种情况能被接受的概率取决于：

$x、T、E(0)$

满足以下公式：
\begin{cases}
x<E(0)：P( E(1)=x ) = 1\
​x>E(0)：P( E(1)=x ) = e^\frac{E(0)-x}{KT}
\end{cases}
看不懂没关系，只需要知道以下结论：

• 可以证明的是，当T趋于无穷大时，最终状态达到全局最小值的概率趋于1。
• 根据公式第二部分很容易看出来，当$x>E(0)$时，仍有概率到达x点。（即达到极小值点后仍有概率跳出极小值，这为得到全局最小值提供了可能，这一点很重要。）
• 在高温情况下（即T较大时），达到高能量状态的概率更大，当温度越来越低时，达到高能量状态的概率将变低。
• 能量的变化量应与温度T有关，很容易理解，当温度降低后发生的能量变化也应该变小。

算法的应用

看了上面关于模型的介绍，你可能还不懂这个算法到达有什么用。我们来简单阐述一下它的应用。当求非线性规划问题的全局最小值时，我们可以使用模拟退火算法进行求解。求解的过程为以下四步：

• 首先任意给一个初始可行解，同时设置一个较大的参数值（相当于模型中的kT）。
• 随机的给出当前点附近的一些可行点（即能量发生变化），它的寻找范围与kT有关。
• 计算可行点对应的值，根据模拟退火模型，比较该点的值与当前值的大小，并通过概率公式判断是否达到该点。
• 不断循环直到达到精度要求。

则最终到达的点有一定概率在最优解附近，再通过贪心的方法即可找到全局最优解。所以这是一个随机算法。

例题

​ 题目来源：BZOJ 3680

题目描述

gty又虐了一场比赛，被虐的蒟蒻们决定吊打gty。gty见大势不好机智的分出了n个分身，但还是被人多势众的蒟蒻抓住了。蒟蒻们将n个gty吊在n根绳子上，每根绳子穿过天台的一个洞。这n根绳子有一个公共的绳结x。吊好gty后蒟蒻们发现由于每个gty重力不同，绳结x在移动。蒟蒻wangxz脑洞大开的决定计算出x最后停留处的坐标，由于他太弱了决定向你求助。
不计摩擦，不计能量损失，由于gty足够矮所以不会掉到地上。

输入数据

输入第一行为一个正整数n(1<=n<=10000)，表示gty的数目。
接下来n行,每行三个整数xi，yi，wi，表示第i个gty的横坐标，纵坐标和重力。
对于20%的数据，gty排列成一条直线。
对于50%的数据，1<=n<=1000。
对于100%的数据，1<=n<=10000,-100000<=xi,yi<=100000

输出数据

输出1行两个浮点数（保留到小数点后3位），表示最终x的横、纵坐标。

样例输入

3

0 0 1

0 2 1

1 1 1

样例输出

0.577 1.000

题目解析

本题用到了数学知识：广义费马点。

简单的说一下道理：因为重物悬挂的位置有尽可能下降的趋势，所以最终稳定的点应为重力势能和的最小点，同时由于绳子长度一定，则重力势能最小对应下拉绳子长度最长，对应剩余绳长最短，即求到所有点的加权距离和最短的点，其中权值为重力。即求广义费马点。

代码

#include<bits/stdc++.h>
using namespace std;
const int maxn = 1e4+5;
int n;
double x[maxn],y[maxn],w[maxn];
double xs=0,ys=0,dis;
double pow2(double x){return x*x;}
double d(double x1,double y1,double x2,double y2){
    return sqrt(pow2(x1-x2)+pow2(y1-y2));
}
double calc(double nx,double ny){
    double res=0;
    for(int i=0;i<n;i++)
        res+=w[i]*d(nx,ny,x[i],y[i]);
    if(res<dis)
        dis=res,xs=nx,ys=ny;
    return res;
}
double Rand(){
    return (double)(rand()%10000)/10000.0;
}
void work(){
    double dE,T=1e5;
    double nowx,nowy,nextx,nexty;
    nowx=xs;nowy=ys;
    //模拟退火
    while(T>0.001){
        nextx=nowx+T*(Rand()*2-1);
        nexty=nowy+T*(Rand()*2-1);
        dE=calc(nowx,nowy)-calc(nextx,nexty);
        if(dE>0 || exp(dE/T)>Rand())
            nowx=nextx,nowy=nexty;
        T*=0.97;
    }
    //在结束退火后找到极小值点
    for(int i=0;i<1000;i++){
        nextx=xs+T*(Rand()*2-1);
        nexty=ys+T*(Rand()*2-1);
        calc(nextx,nexty);
    }
}
int main(){
    srand(1594894);
    scanf("%d",&n);
    for(int i=0;i<n;i++){
        scanf("%lf%lf%lf",&x[i],&y[i],&w[i]);
        xs+=x[i],ys+=y[i];
    }
    xs/=n,ys/=n;
    dis = calc(xs,ys);
    work();
    printf("%.3lf %.3lf\n",xs,ys);
    return 0;
}

参考资料

费马点：http://www.matrix67.com/blog/archives/422?replytocom=922225#respond

模拟退火算法：https://zhuanlan.zhihu.com/p/21277465