模拟退火算法解决TSP问题
算法思想
模拟退火算法(Simulate Anneal,SA)是一种通用概率演算法,用来在一个大的搜寻空间内找寻命题的最优解
模拟退火算法来源于固体退火原理,将固体加温至充分高,再让其徐徐冷却,加温时,固体内部粒子随温升变为无序状,内能增大,而徐徐冷却时粒子渐趋有序,在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小。根据Metropolis准则,粒子在温度T时趋于平衡的概率为e(-ΔE/(kT)),其中E为温度T时的内能,ΔE为其改变量,k为Boltzmann常数。用固体退火模拟组合优化问题,将内能E模拟为目标函数值f,温度T演化成控制参数t,即得到解组合优化问题的模拟退火算法:由初始解i和控制参数初值t开始,对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代,并逐步衰减t值,算法终止时的当前解即为所得近似最优解,这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制,包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。
设计思路
-
初始化温度T,初始解状态S,每个温度t下的迭代次数L;
-
当k = 1,2,……,L时,进行3~6;
-
对当前解进行变换得到新解S’(例如对某些解中的元素进行互换,置换);
-
计算增量Δt′=C(S′)-C(S),其中C(S)为评价函数;
-
若Δt′<0则接受S′作为新的当前解,否则以概率exp(-Δt′/(KT))接受S′作为新的当前解(k为玻尔兹曼常数,数值为:K=1.3806505(24) × 10^-23 J/K);
-
如果满足终止条件则输出当前解作为最优解,结束程序;
-
减小T,转到第2步,直到T小于初始设定的阈值。
程序代码(TSP问题)
第一种
#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include<time.h>
#include<math.h>
#define N 109
#define big 9999999
using namespace std;
double nowDist(double d[N][N],int temp[N],int m){
//计算当前temp路线组合距离
double dist=0;
for(int i=1;i<m;i++)dist+=d[temp[i]][temp[i+1]];
dist+=d[temp[m]][temp[1]];
return dist;
}
int main(){
double d[N][N];//距离矩阵
int temp[N];//临时地点排列顺序表
double dist=0;//最小距离
int n,m;//n为输入行数,m为地点个数
int res[N];
double mostdist=big;
//数据输入
freopen("input.txt","r",stdin);
cin>>m>>n;
while(n-->0){
int i,j;
double k;
cin>>i>>j>>k;
d[i][j]=d[j][i]=k;
}
for(int i=1;i<=m;i++){
d[i][i]=big;
temp[i]=i;
if(i!=m) dist+=d[i][i+1];
}
dist+=d[m][1];
cout<<"初始距离: "<<dist<<endl;
int count=40;
while(count--){//取count次迭代的最小结果
double T=1;//初始温度
double a=0.999;//退火率
double low=1e-30;//最低温度
long ct=0;//迭代步数
long ctj=0;//有效迭代步数
srand(time(0));//随机数种子
long all=300000;//最大迭代次数
while(T>low){
int t1=rand()%m+1,t2=rand()%m+1;//随机获得两个位置
if(t1!=t2){
swap(temp[t1],temp[t2]);//交换两位置
double ndist=nowDist(d,temp,m);//交换后的路线总距离
double df=(ndist-dist);
if(df<0){//当前解更优
dist=ndist;
ctj++;
}else if(exp(-df/T)>(rand()%100)/100.0){//大于
dist=ndist;
ctj++;
}else{
swap(temp[t1],temp[t2]);
}
T*=a;//降温
}
ct++;
if(ct>all)break;
}
cout<<"当前最短距离:"<<dist<<endl;
if(dist<mostdist){
mostdist=dist;
for(int i=1;i<=m;i++)res[i]=temp[i];
}
}
cout<<"最短距离:"<<mostdist<<endl<<"路线:";
for(int i=1;i<=m;i++){
cout<<res[i]<<" ";
}
return 0;
}
第二种
#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include<time.h>
#include<math.h>
#include<fstream>
#include<algorithm>
#include<memory.h>
using namespace std;
const int num = 1000;//city number
const int width = 100;
const int height = 100;
typedef struct node {
int x;
int y;
}city;
city citys[num];//citys
double dic[num][num];//distance from two citys;
bool visit[num];//visited
int N;//real citys
int seq[num];//最优路径序列
double answer;//最优路径长度
const int tempterature = 1000;//初始温度
const double u = 0.998;//成功降温因子
const double v = 0.999;//失败降温因子
int k = 100;//对每个温度迭代次数
void init() {//set N&&x-y
N = 51;
citys[0].x = 37; citys[0].y = 52;
citys[1].x = 49; citys[1].y = 49;
citys[2].x = 52; citys[2].y = 64;
citys[3].x = 20; citys[3].y = 26;
citys[4].x = 40; citys[4].y = 30;
citys[5].x = 21; citys[5].y = 47;
citys[6].x = 17; citys[6].y = 63;
citys[7].x = 31; citys[7].y = 62;
citys[8].x = 52; citys[8].y = 33;
citys[9].x = 51; citys[9].y = 21;
citys[10].x = 42; citys[10].y = 41;
citys[11].x = 31; citys[11].y = 32;
citys[12].x = 5; citys[12].y = 25;
citys[13].x = 12; citys[13].y = 42;
citys[14].x = 36; citys[14].y = 16;
citys[15].x = 52; citys[15].y = 41;
citys[16].x = 27; citys[16].y = 23;
citys[17].x = 17; citys[17].y = 33;
citys[18].x = 13; citys[18].y = 13;
citys[19].x = 57; citys[19].y = 58;
citys[20].x = 62; citys[20].y = 42;
citys[21].x = 42; citys[21].y = 57;
citys[22].x = 16; citys[22].y = 57;
citys[23].x = 8; citys[23].y = 52;
citys[24].x = 7; citys[24].y = 38;
citys[25].x = 27; citys[25].y = 68;
citys[26].x = 30; citys[26].y = 48;
citys[27].x = 43; citys[27].y = 67;
citys[28].x = 58; citys[28].y = 48;
citys[29].x = 58; citys[29].y = 27;
citys[30].x = 37; citys[30].y = 69;
citys[31].x = 38; citys[31].y = 46;
citys[32].x = 46; citys[32].y = 10;
citys[33].x = 61; citys[33].y = 33;
citys[34].x = 62; citys[34].y = 63;
citys[35].x = 63; citys[35].y = 69;
citys[36].x = 32; citys[36].y = 22;
citys[37].x = 45; citys[37].y = 35;
citys[38].x = 59; citys[38].y = 15;
citys[39].x = 5; citys[39].y = 6;
citys[40].x = 10; citys[40].y = 17;
citys[41].x = 21; citys[41].y = 10;
citys[42].x = 5; citys[42].y = 64;
citys[43].x = 30; citys[43].y = 15;
citys[44].x = 39; citys[44].y = 10;
citys[45].x = 32; citys[45].y = 39;
citys[46].x = 25; citys[46].y = 32;
citys[47].x = 25; citys[47].y = 55;
citys[48].x = 48; citys[48].y = 28;
citys[49].x = 56; citys[49].y = 37;
citys[50].x = 30; citys[50].y = 40;
}
void set_dic() {//set distance
for (int i = 0; i<N; ++i) {
for (int j = 0; j<N; ++j) {
dic[i][j] = sqrt(pow(citys[i].x - citys[j].x, 2) + pow(citys[i].y - citys[j].y, 2));
}
}
}
double dic_two_point(city a, city b) {
return sqrt(pow(a.x - b.x, 2) + pow(a.y - b.y, 2));
}
double count_energy(int* conf) {
double temp = 0;
for (int i = 1; i<N; ++i) {
temp += dic_two_point(citys[conf[i]], citys[conf[i - 1]]);
}
temp += dic_two_point(citys[conf[0]], citys[conf[N - 1]]);
return temp;
}
bool metro(double f1, double f2, double t) {
if (f2 < f1)
return true;
//else
// return false;
double p = exp(-(f2 - f1) / t);
int bignum = 1e9;
if (rand() % bignum<p*bignum)
return true;
return false;
}
void generate(int* s) {//随机产生一组新解
bool v[num];
memset(v, false, sizeof(v));
for (int i = 0; i<N; ++i) {
s[i] = rand() % N;
while (v[s[i]]) {
s[i] = rand() % N;
}
v[s[i]] = true;
}
}
void generate1(int* s) {//随机交换序列中的一组城市顺序
int ti = rand() % N;
int tj = ti;
while (ti == tj)
tj = rand() % N;
for (int i = 0; i<N; ++i)
s[i] = seq[i];
swap(s[ti], s[tj]);
}
void generate2(int* s) {//随机交换序列中的两组城市顺序
int ti = rand() % N;
int tj = ti;
int tk = ti;
while (ti == tj)
tj = rand() % N;
while (ti == tj || tj == tk || ti == tk)
tk = rand() % N;
for (int i = 0; i<N; ++i)
s[i] = seq[i];
swap(s[ti], s[tj]);
swap(s[tk], s[tj]);
}
void generate3(int* s) {//随机选序列中的三个城市互相交换顺序
int ti = rand() % N;
int tj = ti;
int tm = rand() % N;
int tn = ti;
while (ti == tj)
tj = rand() % N;
while (tm == tn)
tn = rand() % N;
for (int i = 0; i<N; ++i)
s[i] = seq[i];
swap(s[ti], s[tj]);
swap(s[tm], s[tn]);
}
void generate0(int* s) {//以上三种交换方式等概率选择
int temp = rand() % 3;
if (temp == 0)
generate1(s);
else if (temp == 1)
generate2(s);
else if (temp == 2)
generate3(s);
}
void moni() {
double t = tempterature;
int seq_t[num];
for (int i = 0; i<N; ++i) {//初始化当前序列
seq[i] = seq_t[i] = i;
}
double new_energy = 1, old_energy = 0;
while (t>1e-9&&fabs(new_energy - old_energy)>1e-9) {//温度作为控制变量
int t_k = k;
int seq_tt[num];
while (t_k--&&fabs(new_energy - old_energy)>1e-9) {//迭代次数作为控制变量
generate1(seq_tt);
new_energy = count_energy(seq_tt);//new
old_energy = count_energy(seq_t);//old
if (metro(old_energy, new_energy, t))
for (int i = 0; i < N; ++i)
seq_t[i] = seq_tt[i];
}
new_energy = count_energy(seq_t);//new
old_energy = answer;//old
if (metro(old_energy, new_energy, t)) {
for (int i = 0; i < N; ++i)
seq[i] = seq_t[i];
answer = count_energy(seq);
t *= u;//接受新状态降温因子0.98
}
else
t *= v;//不接受新状态降温因子0.99
}
answer = count_energy(seq);
}
void output() {
cout << "the best road is : \n";
for (int i = 0; i < N; ++i) {
cout << seq[i];
if (i == N - 1)
cout << endl;
else
cout << " -> ";
}
cout << "the length of the road is " << answer << endl;
}
void test() {
ifstream ifile("data.txt");
if (!ifile) {
cout << "open field\n";
return;
}
while (!ifile.eof()) {
int te = 0;
ifile >> te;
ifile >> citys[te - 1].x >> citys[te - 1].y;
}
}
int main() {
srand(time(0));
int t;
while (cin >> t) {//仅作为重启算法开关使用,无意义
init();//使用程序内置数据使用init()函数,
//test();//使用文件读取数据使用test()函数,
set_dic();//计算每个城市之间的距离
moni();//退火
output();//输出
}
return 0;
}
测试例
第一种,通过以下代码生成数据集文件
import random as rd
n=eval(input('城市数量:'))
res=[]
for i in range(1,n):
for j in range(i+1,n+1):
res.append((i,j,rd.random()*100))
f=open('input.txt','w')
f.write(str(n)+' '+str(len(res))+'\n')
for i in res:
f.write(str(i[0])+" "+str(i[1])+" "+str(i[2])+"\n")
f.close()
input('结束')
第二种
1 37 52
2 49 49
3 52 64
4 20 26
5 40 30
6 21 47
7 17 63
8 31 62
9 52 33
10 51 21
11 42 41
12 31 32
13 5 25
14 12 42
15 36 16
16 52 41
17 27 23
18 17 33
19 13 13
20 57 58
21 62 42
22 42 57
23 16 57
24 8 52
25 7 38
26 27 68
27 30 48
28 43 67
29 58 48
30 58 27
31 37 69
32 38 46
33 46 10
34 61 33
35 62 63
36 63 69
37 32 22
38 45 35
39 59 15
40 5 6
41 10 17
42 21 10
43 5 64
44 30 15
45 39 10
46 32 39
47 25 32
48 25 55
49 48 28
50 56 37
51 30 40
运行结果
第一种
初始距离: 587.205
当前最短距离:264.036
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:213.825
当前最短距离:232.311
当前最短距离:278.578
当前最短距离:232.311
当前最短距离:213.825
当前最短距离:216.204
当前最短距离:219.469
当前最短距离:258.466
当前最短距离:219.469
当前最短距离:213.825
当前最短距离:232.311
当前最短距离:278.578
当前最短距离:232.311
当前最短距离:213.825
当前最短距离:216.204
当前最短距离:219.469
当前最短距离:258.466
当前最短距离:219.469
当前最短距离:213.825
当前最短距离:232.311
当前最短距离:278.578
最短距离:213.825
路线:10 2 1 5 8 6 9 7 4 3
第二种
the best road is :
13 -> 24 -> 22 -> 23 -> 50 -> 4 -> 48 -> 9 -> 32 -> 44 -> 38 -> 29 -> 33 -> 20 -> 28 -> 34 -> 35 -> 2 -> 21 -> 27 -> 30 -> 19 -> 1 -> 49 -> 37 -> 8 -> 15 -> 31 -> 10 -> 0 -> 7 -> 25 -> 42 -> 6 -> 47 -> 26 -> 5 -> 45 -> 36 -> 14 -> 43 -> 41 -> 18 -> 39 -> 40 -> 12 -> 16 -> 3 -> 46 -> 11 -> 17
分析
模拟退火算法具有以下优缺点;
1.迭代搜索效率高,并且可以并行化;
2.算法中有一定概率接受比当前解较差的解,因此一定程度上可以跳出局部最优;
3.算法求得的解与初始解状态S无关,因此有一定的鲁棒性;
4.具有渐近收敛性,已在理论上被证明是一种以概率l收敛于全局最优解的全局优化算法。
参考
https://blog.csdn.net/daaikuaichuan/article/details/81381875
https://blog.csdn.net/wordsin/article/details/79915328
https://blog.csdn.net/chenlnehc/article/details/51933802