参考资料:
基本思想
- 初始化:设定初始温度 T T T,初始解状态 S S S,目标函数(能量函数) C ( ⋅ ) C(\cdot) C(⋅),以及每个温度 T T T下的迭代次数 L L L;
- 对于 k = 1 , 2 , . . . , L k=1,2,...,L k=1,2,...,L,循环步骤3-6;
- 在当前解的基础上,产生新解 S ′ S' S′;
- 计算目标函数增量 Δ t ′ = C ( S ′ ) − C ( S ) \Delta t'=C(S')-C(S) Δt′=C(S′)−C(S);
- 基于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′作为当前解;并且,记录当前的最优解;
- 判断是否满足终止条件(温度 T T T降低到设定的阈值以下;或者最优解长时间未得到更新),满足则输出当前最优解;
- 降低温度 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