模拟退火算法C++

模拟退火算法

为什么叫做“退火”,还要从物理学说起

在热力学上,退火(annealing)现象指物体逐渐降温的物理现象,温度愈低,物体的能量状态会低;够低后,液体开始冷凝与结晶,在结晶状态时,系统的能量状态最低。大自然在缓慢降温(亦即,退火)时,可“找到”最低能量状态:结晶。但是,如果过程过急过快,快速降温(亦称「淬炼」,quenching)时,会导致不是最低能态的非晶形。

这里的最低能量状态,也就是我们题目中的最优解。

实现

因为要模拟退火的过程,因此我们先定义一些变量

T:当前温度,由高温到低温,代表算法进行到了什么程度,一般为double类型

ΔT:每次温度的变化率,一般取0.95−0.99,模拟缓慢降温的过程(上一次的温度乘温度变换率即为这一次的温度)

f(x) 当前状态对应的值

上面我们提到,模拟退火会以一定的概率转移到比当前差的解,那么这个概率是多少呢?科学家经过分析,当这个概率为e^(−Δf/T)时最优

那么根据退火的过程,我们不难得到模拟退火的算法流程

  • 枚举温度T
  • 计算出下一步的状态
  • 若下一步的状态比当前状态优或者满足进行转移的条件,进行转移
  • 降温

因为模拟退火算法具有偶然性,因此我们一般需要对一个问题进行多次模拟退火算法

遇到的问题(codevs 1344)

有 N ( <=20 ) 台 PC 放在机房内,现在要求由你选定一台 PC,用共 N-1 条网线从这台机器开始一台接一台地依次连接他们,最后接到哪个以及连接的顺序也是由你选定的,为了节省材料,网线都拉直。求最少需要一次性购买多长的网线。(说白了,就是找出 N 的一个排列 P1 P2 P3 ..PN 然后 P1 -> P2 -> P3 -> … -> PN 找出 |P1P2|+|P2P3|+…+|PN-1PN| 长度的最小值)

输入描述

第一行 N ,下面 N 行,每行分别为机器的坐标(x,y) ( 实数 -100<=x,y<=100 )

输出描述

最小的长度,保留两位小数。

20 个平面上的点,找出一条连接所有的点的线,要求线最短。这种问题被称为最优组合问题。n = 20,传统的方法不再适用,这时候应该采用比较另类的算法,随机化算法。

裸的随机化算法

直接设置循环次数,每次随机打乱20个点的次序,再求长度,每次循环更新最小值。

伪代码

dis = ∞
for i in 999999:
    shuffle_array([p1..p20])
    temp_dis = sum_dis([p1..p20])
    if temp_dis < dis:
        dis = temp_dis
return dis

code

#include<cstdio>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include<algorithm>

const int maxn = 21;
double x[maxn], y[maxn];
double dist[maxn][maxn];
int path[maxn];
int n;
double path_dist(){
    double ans = 0;
    for(int i = 1; i < n; i++) {
        ans += dist[path[i - 1]][path[i]];
    }
    return ans;
}
int main(){
    srand(19260817U);                            // 使用确定的种子初始化随机函数是不错的选择 
    scanf("%d", &n);
    for(int i = 0; i < n; i++) scanf("%lf%lf", x + i, y + i);
    for(int i = 0; i < n; i++) for(int j = i + 1; j < n; j++) dist[i][j] = dist[j][i] = hypot(x[i] - x[j], y[i] - y[j]);

    for(int i = 0; i < n; i++) path[i] = i;        // 获取初始排列 
    double ans = path_dist();                    // 初始答案 
    int T = 30000000 / n;                         // 单次计算的复杂度是O(n),这里的30000000是试出来的 
    while(T--){
        std::random_shuffle(path, path + n);    // 随机打乱排列 
        ans = std::min(ans, path_dist());        // 更新最小值 
    }
    printf("%.2lf", ans);
}

用这种方法过掉了系统的样例以及5组后台数据。但这种完全随机方法的问题出现在,假如某一次计算尝试出了一个比较好的路径,那么最优的路径很可能可以在原基础上作一两次改动就可以得到,这时候完全打乱整个序列不是一个很好的选择。

爬山法

基于贪心策略的爬山法,每次只顾眼前。根据原序列生成一个新的序列,然后交换新序列的任意两个数。假如说新生成的序列更优,则使用新序列继续计算,否则序列不变。
《模拟退火算法C++》
如图,以爬山举例,假如使用这种方法,我们的起点在C,则我们的算法会带着我们走向A(局部最优解),但由于只顾眼前,发现走向E不会使当前变得更优秀,就不会去尝试。而到达不了B(全局最优解)。
当然,具体的策略还可以通过不断随机初始化序列的排列,再使用这种爬山法求取各处的局部最优解,也能达到全局最优,这里不再赘述。

模拟退火算法

“模拟退火算法以一定的概率来接受一个比当前解要差的解,因此有可能会跳出这个局部的最优解,达到全局的最优解。以图1为例,模拟退火算法在搜索到局部最优解A后,会以一定的概率接受到E的移动。也许经过几次这样的不是局部最优的移动后会到达D点,于是就跳出了局部最大值A。”

“我们要找到山脉的最高峰,但是我(计算机)只能看到我的脚下哪边是上升的,哪边是下降的,看不到远处是否上升。每次移动,我们随机选择一个方向。如果这个方向是上升的的(更优),那么就决定往那个方向走;如果这个方向是下降的(更差),那么“随机地接受”这个方向,接受就走,不接受就再随机一次——这个随机是关键,要考虑很多因素。比如,一个陡的下坡的接受率要比一个缓的下坡要小(因为陡的下坡后是答案的概率小);同样的下降坡度,接受的概率随时间降低(逐渐降低才能趋向稳定)。”

既然这个随机的接受这么重要,那么我们就将它写为一个函数:

bool accept(double delta, double temper){
    if(delta <= 0) return true;
    return rand() <= exp((-delta) / temper) * RAND_MAX;
} 

delta 是新序列长度与原序列长度的差,temper是设置的温度,逐渐降低。
如果delta <=0 说明新序列的长度更短,则返回true接受。
如果delta >0 也就是局部较糟糕的选择,我们随机的返回true。
这里把return rand() <= exp((-delta) / temper) * RAND_MAX;
恢复成移项前的公式
return (double) rand()/RAND_MAX <= exp((-delta) / temper)
在右边, temper是正数,delta是正数(因为delta是负数的已经return出去了),因此exp中间的参数是负数,我们知道指数函数在参数是负数的时候返回(0,1), 我们在左边随机一个实数,如果它比概率小,就接受,否则就不接受。

模拟退火的随机接受的概率逐渐降低是怎么体现的?
主函数中,使temper这个变量逐渐变小。使得 delta / temper 逐渐变大,使得(-delta) / temper 逐渐变小,使得 exp((-delta) / temper) 逐渐变小,使得 左边 <= 右边 的几率逐渐变小。
坡度越陡,接受的概率越小?
坡度越陡,也就是delta越大,使得 delta / temper 逐渐变大,同理使得左边 <= 右边 的几率逐渐变小。

完整代码

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <cmath>
#include <cmath>
#include <algorithm>
using namespace std;
const int maxn = 21;
double x[maxn],y[maxn];
double dist[maxn][maxn];
int n;
const double INF = 0x3f3f3f3f;

/*
 * 模拟退火算法,核心是, 随机的接受,接受的概率随时间降低。(逐渐降低才能趋向于稳定)
 * */
struct Path{
    int path[maxn];
    Path(){
        for(int i=0;i<n;i++) path[i] = i;
    }
    Path(const Path &p){
        memcpy(path, p.path,sizeof path);
        swap(path[rand() % n], path[rand() % n]);

    }
    double dist(){
        double ans = 0;
        for(int i=1;i <n;i++){
            ans += ::dist[path[ i -1]][path[i]];
        }
        return ans;
    }
};

bool accept(double delta, double temper){
    if(delta <= 0) return true;
    return rand() <= exp( (-delta) / temper) * RAND_MAX;
    /*
     * 这里的移项之后观察
     * return (double) rand()/RAND_MAX <= exp((-delta) / temper)
     *
     * 在右边, temper是正数,delta是正数(因为delta是负数的已经return出去了),因此exp中间的参数是负数,我们知道
     * 指数函数在参数是负数的时候返回(0,1), 我们在左边随机一个实数,如果它比概率小,就接受,否则就不接受。
     *
     * 模拟退火的随机接受的概率逐渐降低是怎么体现的?
     *
     * 主函数中,使temper这个变量逐渐变小。
     * 使得 delta / temper 逐渐变大,使得(-delta) / temper 逐渐变小,
     * 使得 exp((-delta) / temper) 逐渐变小,使得 左边 <= 右边 的几率逐渐变小。
     *
     * */
}

double solve(){
    const double max_temper = 10000;
    const double dec = 0.999;
    double temp = max_temper;
    Path p;
    while(temp > 0.01){
        Path p2(p);
        if (accept(p2.dist() - p.dist(),temp)) p = p2;
        temp *= dec;
    }
    return p.dist();

}

int main(){
    srand(19260817U);
    cin >> n;
    for(int i = 0; i < n; i++) {
        scanf("%lf%lf", x + i, y + i);
    }
    for(int i = 0; i < n; i++){
        dist[i][i] = 0;
        for(int j = i + 1; j < n; j++){
            dist[i][j] = dist[j][i] = hypot(x[i] - x[j], y[i] - y[j]);
        }
    }
    double ans = INF;
    int T = 155;
    while(T--){
        ans = min(ans, solve());
    }
    printf("%.2lf", ans);
}


/*
 * 随机化优化算法, 感谢博主https://www.cnblogs.com/vactor/p/8824328.html 的介绍,与思路分享。
 *
 *  solution1. 随机化算法
 * */

参考资料
大白话介绍模拟退火算法
模拟退火算法c++
模拟退火算法

点赞
  1. 苏叶biu说道:

    谢谢老柴,看懂了