模拟退火算法(粗略了解)

参考资料:

基本思想

  1. 初始化:设定初始温度 T T T,初始解状态 S S S,目标函数(能量函数) C ( ⋅ ) C(\cdot) C(),以及每个温度 T T T下的迭代次数 L L L
  2. 对于 k = 1 , 2 , . . . , L k=1,2,...,L k=1,2,...,L,循环步骤3-6;
  3. 在当前解的基础上,产生新解 S ′ S' S;
  4. 计算目标函数增量 Δ t ′ = C ( S ′ ) − C ( S ) \Delta t'=C(S')-C(S) Δt=C(S)C(S);
  5. 基于Metropolis准则判断新解是否被接受:如果 Δ t ′ < 0 \Delta t'<0 Δt<0,直接接受 S ′ S' S作为当前解;否则,以概率 e x p ( − Δ t ′ / T ) exp(-\Delta t'/T) exp(Δt/T)接受新解 S ′ S' S作为当前解;并且,记录当前的最优解;
  6. 判断是否满足终止条件(温度 T T T降低到设定的阈值以下;或者最优解长时间未得到更新),满足则输出当前最优解;
  7. 降低温度 T T T(退火操作): T = δ ⋅ T , 0 < δ < 1 T = \delta\cdot T, 0 < \delta < 1 T=δT,0<δ<1,转回第2步。

Tips

  • 新解的产生:通常选择有当前解通过简单变换即可产生新解的方法,便于后续目标函数计算,减少算法耗时;
  • 目标函数计算:结果变化由新解相对与当前解的变换部分导致,因此最好采取增量式计算新解与当前解的目标函数差值(如果目标函数十分复杂,可以极大的减少计算量)。

C++示例代码(TSP问题)

输入数据 文件TSP.data

6
20 80
16 84
23 66
62 90
11 9
35 28

代码(目标函数比较简单,因此也并没有增量式计算目标函数差值)

// 主体来源:http://www.imooc.com/article/details/id/255581
// 修改了SA算法实现部分
#include <iostream>
#include <string>
#include <cstdlib>
#include <algorithm>
#include <cstdio>
#include <ctime>
#include <cmath>

#define N     30      //城市数量
#define T     3000    //初始温度
#define EPS   1e-7    //终止温度,用于退出条件
#define DELTA 0.98    //温度衰减率,用于退出条件

#define ILOOP 1000   //同一温度T下循环次数,即步数L
#define MAXSTAYCOUNTER 150 //最优解未更新上限,用于判断退出条件

using namespace std;

//定义路线结构体
struct Path
{
    int citys[N];
    double len;
};

//定义城市点坐标
struct Point
{
    double x, y;
};

Path bestPath;        //记录最优路径
Point p[N];       //每个城市的坐标
double w[N][N];   //两两城市之间路径长度
int nCase;        //测试次数

// 计算两点距离
double dist(Point A, Point B)
{
    return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}

// 计算载入城市之间距离
void GetDist(Point p[], int n)
{
    for(int i = 0; i < n; i++){
        for(int j = i + 1; j < n; j++){
            w[i][j] = dist(p[i], p[j]);
            w[j][i] = w[i][j];
        }
    }
}

// 从文件中读取数据
void Input(Point p[], int &n)
{
    scanf("%d", &n);
    for(int i = 0; i < n; i++)
        scanf("%lf %lf", &p[i].x, &p[i].y);
}

// 初始化解,按城市顺序
void Init(int n)
{
    nCase = 0;
    bestPath.len = 0;
    for(int i = 0; i < n; i++)
    {
        bestPath.citys[i] = i;
        if(i != n - 1)
        {
            printf("%d--->", i);
            bestPath.len += w[i][i + 1];
            cout << "(" << bestPath.len << ")";
        }
        else
            printf("%d\n", i);
    }
    cout << "\nInit path length is : " << bestPath.len << endl;
    cout << "-----------------------------------" << endl << endl;
}

// 输出解t
void Print(Path t, int n)
{
    printf("Path is : ");
    for(int i = 0; i < n; i++)
    {
        if(i != n - 1)
            printf("%d-->", t.citys[i]);
        else
            printf("%d\n", t.citys[i]);
    }
    cout << "\nThe path length is : " << t.len << endl;
    cout << "-----------------------------------" << endl << endl;
}

// 产生新解
Path GetNext(Path p, int n)
{
    Path ans = p;
    int x = (int)(n * (rand() / (RAND_MAX + 1.0)));
    int y = (int)(n * (rand() / (RAND_MAX + 1.0)));
    while(x == y)
    {
        x = (int)(n * (rand() / (RAND_MAX + 1.0)));
        y = (int)(n * (rand() / (RAND_MAX + 1.0)));
    }
    swap(ans.citys[x], ans.citys[y]);
    ans.len = 0;
    for(int i = 0; i < n - 1; i++)
        ans.len += w[ans.citys[i]][ans.citys[i + 1]];
    // cout << "nCase = " << nCase << endl;
    // Print(ans, n);
    nCase++;
    return ans;
}

// 模拟退火,核心算法,
// 更改自Python版本:https://github.com/guofei9987/scikit-opt/blob/master/sko/SA.py
void SA_Curya(int n){
    double t = T;
    srand((unsigned)(time(NULL)));
    Path curPath = bestPath;
    Path newPath = bestPath;
    int iter_cycle = 0;
    int stay_counter = 0;

    while(true){
        double curBest = bestPath.len;
        for(int i = 0; i < ILOOP; ++i){
            // 产生一个新解
            newPath = GetNext(curPath, n);
            double dE = newPath.len - curPath.len;

            // Metropolis,判断新解是否被接受
            if(dE < 0 || exp(-dE / t) > (rand() / (RAND_MAX + 1.0))){
                curPath = newPath;
                if(newPath.len < bestPath.len){
                    bestPath = newPath;
                    Print(bestPath, n);
                }
            }
        }
        // 退火降温,同时计数器加一
        t *= DELTA;
        iter_cycle++;

        // 记录最优解未变化的迭代次数
        if(abs(bestPath.len-curBest) <= max(1e-9*max(abs(bestPath.len), abs(curBest)), 1e-30))
            stay_counter++;
        else
            stay_counter = 0;

        // 退出条件:
        // 1、温度降低到最小值
        if(t < EPS){
            cout << "Cooled to final temperature!" << endl;
            break;
        }
        // 2、最优解长时间未更新
        if(stay_counter > MAXSTAYCOUNTER){
            cout << "Stay unchanged in " << stay_counter << " iterations!" << endl;
            break;
        }
    }
}

int main(int argc, const char * argv[]) {

    freopen("TSP.data", "r", stdin);
    int n;
    Input(p, n);
    GetDist(p, n);
    Init(n);
    SA_Curya(n);
    Print(bestPath, n);
    printf("Total test times is : %d\n", nCase);
    return 0;
}

输出

0--->(5.65685)1--->(24.9701)2--->(70.7631)3--->(166.481)4--->(197.092)5

Init path length is : 197.092
-----------------------------------

Path is : 0-->1-->2-->3-->5-->4

The path length is : 168.997
-----------------------------------

Path is : 2-->1-->0-->3-->5-->4

The path length is : 166.379
-----------------------------------

Path is : 4-->5-->3-->1-->0-->2

The path length is : 164.599
-----------------------------------

Path is : 4-->5-->3-->2-->0-->1

The path length is : 164.002
-----------------------------------

Path is : 3-->0-->1-->2-->5-->4

The path length is : 138.604
-----------------------------------

Path is : 3-->1-->0-->2-->5-->4

The path length is : 136.825
-----------------------------------

Stay unchanged in 151 iterations!
Path is : 3-->1-->0-->2-->5-->4

The path length is : 136.825
-----------------------------------

Total test times is : 152000
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值