【学习笔记】模拟退火

• 参考：

        模拟退火目的是求一个函数f(x)的最值，这个函数可能不具有任何特殊的单调性质；

          http://www.cnblogs.com/rvalue/p/8678318.html

       http://www.cnblogs.com/flashhu/p/8884132.html

---

• 爬山算法： 

        先随机选取一个值作为初始的$x$，然后每次随机一个$Delta$,$x' = x + Delta x$，如果$f(x')$比$f(x)$更优，那么直接转移到x'，否则不转移;

       很容易看出这样的算法的缺陷在于会陷入局部最优解，想像一个波浪状的函数，很有可能在波峰周围的值都比它小但在远处还存在另一个更高的波峰，需要改进使得算法更加可靠

---

• 理论：

        模拟退火参考了晶体的退火理论(我也不是很专业的说，只能大致表明其意思)，在晶体冷却的过程中，最后的稳定状态宏观是内能最低的，设$Delta E$为内能变化，$T$为当前温度并在不断变低，$k$为一个正常数，则发生状态转移的概率为：$$P(Delta E) = e^{frac{Delta E}{kT} }$$

       以寻求最大值为例，把$Delta E$看成$f(x') - f(x)$ , 带入上述公式可以发现，当$x'$比$x$更优时$f(x') - f(x)>0$，$P(Delta E)  > 1$所以一定接受$x'$ ；当$x$不比$x'$更优时，根据上述公式$P(Delta E) lt 1$，并且随着温度减小，对于相同的$Delta E$ , $P(Delta E)$在减小。用来修改爬山算法相当于：我们并不以是否更优作为接受几个新的值的标准，而是给出了一个限度，并且这个限度随着温度$T$的减小不断收紧直到确定了答案。

---

• 例子及参数调整：

       用算法去实现需要设计参数:$T0，x0$表示最初温度和最初值，$T$表示当前温度，$T_{k}$表示稳定温度，$d$表示每次的温度变化率，$A$表示每次生成范围所需的常数，可以如下描述模拟退火：

      

T = T0 , x = x0;
while(T > Tk){
    x' = rand(A*T,-A*T);
    delta E = f(x') - f(x);
    P0 = e^(delta E / (k * T) ) ;
    P = rand(0,1); 
    if(P < P0)  x = x' ; 
    T = T  *  d;  
}

　　调参技巧：

     ① A*T是为了让每次的范围和温度成正比最后趋于稳定，一般可以取A=1，这个参数也可以调的；

     ②x0可选取一个可能的最优值，T一般较大，Tk一般取一个较小的正数，d一般选取一个接近1的正数0.980往上走之类的；

     ③T0,Tk和d可以决定调整的次数，即精度和准确度，A可以决定跳出局部解的能力，这些值需要比较精细地测试，次数太少，容易跑不出最优解而陷入局部解，次数太多又容易乱跳；但是在实践中，可以记录每次退火的所有x的最优值有效避免次数过多的问题，所以可以尽量开大再不断调整；

       ④可以用while(clock()<CLOCKS_PER_SEC*0.9) 卡卡时间多次退火取最优 ， 活用rand()函数和RAND_MAX ; 

     ⑤最大值和最小值在符号处理上有点区别需要注意一下；

      给一个例子：

      bzoj3680 吊打XXX

      问题即是使得选定点到所有点的带权距离和最小，可以退火实现

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<queue>
#include<cmath>
#include<vector>
#include<stack>
#include<map>
#include<set>
#include<ctime>
#define Run(i,l,r) for(int i=l;i<=r;i++)
#define Don(i,l,r) for(int i=l;i>=r;i--)
#define ll long long
#define ld double
#define inf 1e18 
#define il inline 
#define Max 10000
#define rd (rand()*2-Max)
using namespace std;
const int N=1e5;
int n,px[N],py[N],pw[N];
ld tk=1e-14,t=2000,d=0.980,ansx,ansy,nowx,nowy,ans,now;
il ld cal(ld x,ld y){
    ld re=0;
    for(int i=1;i<=n;i++){
        re+=1.0* pw[i] * sqrt( (px[i]-x)*(px[i]-x) + (py[i]-y)*(py[i]-y) );
    }
    return re; 
}
int main(){
//  freopen("bzoj3680.in","r",stdin);
//  freopen("bzoj3680.out","w",stdout);
    srand(time(NULL));
    scanf("%d",&n);
    Run(i,1,n)scanf("%d%d%d",&px[i],&py[i],&pw[i]),ansx+=px[i],ansy+=py[i];
    ansx/=n,ansy/=n;
    ans=cal(ansx,ansy);
    while(t>tk){
        nowx = ansx + rd*t;//(RD*2-1)*t; 
        nowy = ansy + rd*t;//(RD*2-1)*t;
    //  nowx = ansx + ((ld)rd/Max*1000)*t; 
    //  nowy = ansy + ((ld)rd/Max*1000)*t;
        now = cal(nowx,nowy);
        ld D = now - ans;
        if(D<0)ansx=nowx,ansy=nowy,ans=now;
        else if(exp(-D/t)*Max>rand())ansx=nowx,ansy=nowy,ans=now; 
        t=t*d;
    }
    printf("%.3lf %.3lf
",ansx,ansy);
    return 0;
}//by tkys_Austin;
 

      另外对于排列相关的问题例如TSP也是适用的

    bzoj2428[HAOI2006]均分数据

    然而我用的随机贪心

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<queue>
#include<cmath>
#include<vector>
#include<stack>
#include<map>
#include<set>
#define Run(i,l,r) for(int i=l;i<=r;i++)
#define Don(i,l,r) for(int i=l;i>=r;i--)
#define ll long long
#define ld long double
#define inf 0x3f3f3f3f
using namespace std;
const int N=21;
int n,m,f[N][N],ans,a[N],now,tmp[N],s[N];
/*
int dp(){
    memset(f,0x3f,sizeof(f));
    for(int i=1;i<=n;i++)s[i]=a[i]+s[i-1];
    f[0][0]=0;
    for(int i=1;i<=n;i++)
    for(int j=1;j<=m;j++){
        for(int k=0;k<j;k++)f[i][j]=min(f[i][j],f[k][j-1]+(s[i]-s[k])*(s[i]-s[k]));
    }
    return f[n][m];
}*/
int cal(){
    int re=0;
    for(int i=1;i<=m;i++)tmp[i]=0;
    for(int i=1;i<=n;i++){
        int p = min_element(tmp+1,tmp+m+1) - tmp;
        tmp[p]+=a[i];
    }
    for(int i=1;i<=m;i++)re+=tmp[i]*tmp[i];
    return re;
}
int main(){
//    freopen("bzoj2428.in","r",stdin);
//    freopen("bzoj2428.out","w",stdout);
    scanf("%d%d",&n,&m);
    ld sum=0;
    for(int i=1;i<=n;i++)scanf("%d",&a[i]),sum+=a[i];
    ans = cal();
    for(int i=1;i<=500000;i++){
        random_shuffle(a+1,a+n+1);
        ans = min(ans , cal());
    }
    sum=sum*sum/m;
    printf("%.2Lf
",sqrt((ans-sum)/m));
    return 0;
}//by tkys_Austin;