模拟退火

一、模拟退火

模拟退火算法用于计算运算量大或随机概率较高的多峰函数最值问题，在多次退火下，正确的概率表现的还是非常出色的。

1.1 \(\quad\) 劣解的接受与 Metropolis 准则

爬山算法只能应用于单峰函数，因为它每次只在附近只寻找更加优秀的解。放在多峰函数下容易陷入局部最大值，而无法找到全局最大值。

不同于爬山算法，在模拟退火算法中，我们在当前位置的一定范围内随机一个位置进行决策。如果这个决策比现在的决策更加优秀，我们无条件地接受；如果这个决策不如当前的决策我们以某种概率接受这个劣解。

具体地，我们像冶金工业退火一样，一开始，我们有一个温度 \(T\)，表示当前的活跃性。这个温度随着随机次数的增加而降低，当最终小于某一个温度 \(t_0\) 时就终止退火。

我们设接受一个比当前解劣 \(\Delta E\) 的劣解的概率为 \(P(\Delta E)\)。根据 Metropolis 准则，我们划定这个概率，并让其与当前温度有关。即表示：随机次数较小时，我们有更大概率接受劣解；随机次数过多时，我们有较小的概率接受劣解。这样既能保证向着最大值的方向查找，又能避免陷入局部最大值。

Metropolis 准则 \(\quad\) 当前温度为 \(T\) 时，对于新状态 \(S'\) 与当前最优状态 \(S\) 的关于最优值的差为 \(\Delta E\ge 0\)，则发生状态转移（接受新状态）的概率为： \[ P(\Delta E)= \begin{cases} 1,&S'\text{比} S \text{更优,} \\ e^{\frac{-\Delta E}{T}}, & \text{otherwise.} \end{cases} \]

1.2 \(\quad\) 算法的实现

一开始，我们有三个参数：

  • 初始温度 \(T_0\)；
  • 降温系数 \(d\)，即每经过依次随机，温度变为 \(T_0d\)，\(d\) 一般为趋近于 \(1\) 的小数；
  • 终止温度 \(T_k\)。

对于每次随机，我们在当前最优解的附近随机新状态（随机范围也与当前温度有关），经计算后依靠 Metropolis 准则决定是否接受。

在依靠 Metropolis 准则决策时需要注意：

  • 在 C++ 语言中，我们可以用 cmath 库中的 exp 函数计算以 \(e\) 为底的指数函数。具体地，可以用 exp(-E/T) 表示 \(e^{\frac{-\Delta E}{T}}\)；
  • 注意计算接受劣解的概率时 \(e\) 的指数是负数；
  • 根据函数图像可知，对于函数 \(f(x)=e^x(x<0)\) 的值域是 \((0,1)\)。所以我们可以将 \(e^{\frac{-\Delta E}{T}}\) 与一个在 \((0,1)\) 范围内的随机数比较。需要注意的是，因为随机次数越多，\(T\) 越小，\(\frac{-\Delta E}{T}\) 也越小，\(f(x)\) 的函数值也越趋近于 \(0\)，接受的概率也应该越低。所以，无论求最大值还是最小值，都只能依靠 exp(-E/T) >= randFromRange(0,1) 或其他单调判断来决策，注意不等号方向和随机值范围，否则有概率退化成爬山算法。其中 randFromRange 函数是在 \((l,r)\) 范围内生成随机数。

一般利用模拟退火解决问题，通常有如下设置：

  • \(T_0\) 一般取 \([2000,3000]\) 中的数；
  • \(d\) 通常取 \(0.999\)，这个值取决于不同题目；
  • \(T_k\) 通常取 \([10^{-16},10^{-5}]\) 之间；
  • 通常情况下会进行多轮退火。每轮结束后不必要清空最优值，只需要重置温度等参数即可。

二、例题练习

2.1 \(\quad\) 分金币

题目来源：TJOI2010，洛谷P3878

例题 \(\quad\) 现在有 \(n\) 枚金币，第 \(i\) 枚金币的价值为 \(v_i\)。现在要把它们分成两部分，要求这两部分金币数目之差不超过 \(1\)，求这样分成的两部分金币的价值之差的最小值。

考虑模拟退火，先把原序列随便分成两个部分，之后随机交换两个部分中的两个数。

核心代码如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
mt19937 rnd(time(0));
void sa(){
    double temp=2333;//初始温度
    while(temp>=1e-15){//终止温度
        int x=rnd()%mid+1;//在两个部分随机一个数交换
        int y=rnd()%(n-mid)+mid+1;
        ll dis=abs(s1-a[x]+a[y]-(s2-a[y]+a[x]));//计算新贡献
        ll D=dis-ans;
        if(D<0){//优解，注意D和上文定义的E是相反数，D是负数表示的是解较优
            ans=dis;
            s1=s1-a[x]+a[y];
            s2=s2-a[y]+a[x];
            swap(a[x],a[y]);
        }
        else if(exp(D/temp)>1.*(rnd()%100000)/100000){//劣解，注意随机范围和不等号方向
            s1=s1-a[x]+a[y];
            s2=s2-a[y]+a[x];
            swap(a[x],a[y]);
        }
        temp*=0.999;//每次温度降低
    }
}
int main(){
    int T;
    cin>>T;
    while(T--){
        cin>>n;
        for(int i=1;i<=n;i++) cin>>a[i];
        mid=(1+n)>>1;
        s1=s2=0;
        for(int i=1;i<=mid;i++) s1+=a[i];
        for(int i=mid+1;i<=n;i++) s2+=a[i];
        ans=abs(s1-s2);
        if(n==1){
            cout<<ans<<endl;
            continue;
        }
        for(int i=1;i<=20;i++) sa();//多进行几轮退火
        cout<<ans<<endl;
    }
    return 0;
}

2.2 \(\quad\) Run Away

题目来源：SP34

例题 \(\quad\) 给定 \(n\) 个点的坐标，在给定范围内找一个点，使得距离所有点的最小值最大。

考虑模拟退火，一开始将这个点的位置设置在范围的正中心。每次在一定范围内随机新点，范围与当前温度有关，见代码。

核心代码如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
void sa(){
    double temp=10000;
    while(temp>=1e-4){
        //随机新点
        double nowx=getrnd(max(0.,nx-temp),min(1.*X,nx+temp));
        double nowy=getrnd(max(0.,ny-temp),min(1.*Y,ny+temp));
        //计算新值
        double dis=calc(nowx,nowy);
        double D=dis-ans;
        //方案较优
        if(D>0){
            ans=dis;
            ansx=nx=nowx;
            ansy=ny=nowy;
        }
        //以一定概率接受劣解，注意不等号方向
        else if(exp(D/temp)>=getrnd(0,1)){
            nx=nowx;
            ny=nowy;
        }
        temp*=0.999;
    }
}
int main(){
    cin>>T;
    while(T--){
        cin>>X>>Y>>n;
        for(int i=1;i<=n;i++) cin>>ax[i]>>ay[i];
        ansx=nx=1.*X/2;
        ansy=ny=1.*Y/2;
        ans=calc(nx,ny);
        for(int i=1;i<=100;i++) sa();//多进行几次退火
        printf("The safest point is (%.1lf, %.1lf).\n",ansx,ansy);
    }
    return 0;
}