模拟退火学习

模拟退火，很奇妙！

原理

现在要求答案 \(f(S)\) 的最值。（假设要求最小值，函数为实值函数。）
设有一个 \(S'\) 是与 \(S\) 接近的状态。不妨令 \(\Delta=f(S')-f(S)\)。那么更新答案让 \(S\) 变成 \(S'\) 的概率为：

\[P(\Delta)=\left\{\begin{matrix}1&\Delta\le0\\e^{-\frac\Delta T}&\Delta> 0\end{matrix}\right. \]

其中温度为 \(T\)，让 \(T\) 逐渐减小，即是 \(T\leftarrow T\times d\)，直到 \(T\le T_k\)。
若要求最大值，将公式中的 \(-\Delta\) 改成 \(\Delta\) 即可。
通常取 \(T=10^9\)，\(T_k=10^{-9}\)，\(d=0.999\)。

代码

随机数

inline double rnd(double l,double r){
    return (double)rand()/RAND_MAX*(r-l)+l;
}
inline int rnd(int l,int r){
    return rand()%(r-l+1)+l;
}

模退多次

void work(){
    srand(time(0));
    for(;(double)clock()/CLOCKS_PER_SEC<limit;)sa(T,Tk,d);
}

接下来是不同的模退类型。
小技巧：\(P(\Delta)=P(e^{-\frac\Delta T}>rand(0,1))\)。其中 \(rand(0,1)\) 为在 \(0,1\) 范围内随机生成的实数。

模退 \(S\) 为实数

double sa(double T,double Tk,double d){
    double pos=rnd(-A,A),val=calc(pos),ans=val;
    for(;T>Tk;T*=d){
        double nwpos=rnd(pos-T,pos+T),nwval=calc(nwpos),delta=nwval-val;
        ans=min(ans,nwval);
        if(exp(-delta/T)>rnd(0,1))pos=nwpos,val=nwval;
    }
    return ans;
}

模退 \(S\) 为向量（点）

double sa(double T,double Tk,double d){
    Point pos=Point(rnd(-A,A),rnd(-A,A));
    double val=calc(pos),ans=val;
    for(;T>Tk;T*=d){
        Point nwpos=Point(rnd(pos.x-T,pos.x+T),rnd(pos.y-T,pos.y+T));
        double nwval=calc(nwpos),delta=nwval-val;
        ans=min(ans,nwval);
        if(exp(-delta/T)>rnd(0,1))pos=nwpos,val=nwval;
    }
    return ans;
}

模退 \(S\) 为排列（排列为 \(1\) 到 \(n\) 的排列）

double sa(double T,double Tk,double d){
    for(int i=1;i<=n;i++)p[i]=i;
    random_shuffle(p+1,p+n+1);
    double val=calc(),ans=val;
    for(;T>Tk;T*=d){
        int pos1=rnd(1,n),pos2=(1,n);
        swap(p[pos1],p[pos2]);
        double nwval=calc(),delta=nwval-val;
        ans=min(ans,nwval);
        if(exp(-delta/T)>rnd(0,1))val=nwval;
        else swap(p[pos1],p[pos2]);//相当于没有操作
    }
    return ans;
}

模退 \(S\) 为集合（集合元素为 \(1\) 到 \(n\) 的整数）

double sa(double T,double Tk,double d){
    for(int i=1;i<=n;i++)v[i]=rnd(0,1);
    double val=calc(),ans=val;
    for(;T>Tk;T*=d){
        int pos=rnd(1,n);
        v[pos]^=1;
        double nwval=calc(),delta=nwval-val;
        ans=min(ans,nwval);
        if(exp(-delta/T)>rnd(0,1))val=nwval;
        else v[pos]^=1;
    }
    return ans;
}