前几天在做孔群加工问题，各种假设到最后就是求解旅行商问题了，因为原本就有matlab代码模板所以当时就改了城市坐标直接用了，发现运行速度惨不忍睹，最后用上了两个队友的电脑一起跑。这次模拟结束后在想用c语言来实现的话应该可以提高不少效率。关于模拟退火和旅行商问题的介绍我就不赘述了，网上各路大神说的都很详细，我下面就把c语言和matlab代码先附上。
c语言：
#ifndef _OPTION_H #define _OPTION_H /* * T0 表示 初始温度 * Tf 表示 结束时的温度 * alpha 表示 温度衰减系数 * T 表示 当前温度 * Markov_length 表示 Markov链长度 */ ; double alpha = 0.99; ; ; double T = T0; #endif // !_OPTION_H
OPTION.h
#ifndef _CITY_H
#define _CITY_H
struct CITY
{
int id;
double x;
double y;
};
#endif
CITY.h
#include <cstdio>
#include <cmath>
#include <windows.h>
#include <ctime>
#include <algorithm>
#include "OPTION.h"
#include "CITY.h"
using namespace std;
typedef unsigned int UINT;
typedef long long LL;
<< );
;
][]; /* 城市之间的距离矩阵 */
]; /* 产生的新解 */
]; /* 当前解 */
]; /* 冷却中得到的最优解 */
double new_dist; /* 新解的路径距离 */
double cur_dist; /* 当前解对应路径距离 */
double best_dist; /* 冷却过程中得到的最优解 */
CITY citys[]; /* 保存城市坐标和编号信息 */
; /* 保存城市数量 */
char *file = NULL; /* 数据文件 */
time_t st, ed; /* 记录开始结束时间 */
void put_help();
double init_rand();
void input();
void cal_dist();
void init();
void cool();
void disp_result();
int main(int argc, char **argv)
{
srand(static_cast<int>(time(NULL)));
) {
put_help();
goto end;
}
/* get the parameters */
file = argv[];
)
{
; i < argc; i += ) {
) {
T0 = atof(argv[i + ]);
T = T0;
}
) {
Tf = atof(argv[i + ]);
}
){
alpha = atof(argv[i + ]);
}
) {
Markov_length = atoi(argv[i + ]);
}
}
}
freopen(file, "r", stdin);
input();
st = clock();
cal_dist();
init();
cool();
ed = clock();
disp_result();
end:
;
}
double init_rand()
{
return rand() / rand_max;
}
void input()
{
double x, y;
while (~scanf("%lf%lf", &x, &y)) {
++city_count;
citys[city_count].id = city_count;
citys[city_count].x = x;
citys[city_count].y = y;
}
puts("Data input success!");
}
void cal_dist()
{
; i <= city_count; ++i) {
; j <= city_count; ++j) {
dist[i][j] = sqrt((citys[i].x - citys[j].x)*(citys[i].x - citys[j].x)
+ (citys[i].y - citys[j].y)*(citys[i].y - citys[j].y));
}
}
}
void init()
{
cur_dist = inf;
best_dist = inf;
; i <= city_count; ++i) {
new_sol[i] = i;
cur_sol[i] = i;
best_sol[i] = i;
}
}
void cool()
{
while (T >= Tf) {
; i <= Markov_length; ++i) {
;
;
|| ind2 == ) {
ind1 = static_cast<int>(ceil(city_count*init_rand()));
ind2 = static_cast<int>(ceil(city_count*init_rand()));
}
swap(new_sol[ind1], new_sol[ind2]);
/* 检查是否满足约束 */
/* 计算距离 */
new_dist = ;
; i <= city_count - ; ++i) {
new_dist = new_dist + dist[new_sol[i]][new_sol[i + ]];
}
new_dist += dist[new_sol[city_count]][new_sol[]];
/* 若新解距离小于当前最优距离,则接受新解；否则以一定概率接受新解 */
if (new_dist < cur_dist) {
cur_dist = new_dist;
; i <= city_count; ++i)
cur_sol[i] = new_sol[i];
if (new_dist < best_dist) {
best_dist = new_dist;
; i <= city_count; ++i)
best_sol[i] = new_sol[i];
}
}
else {
if (init_rand() < exp(-(new_dist - cur_dist) / T)) {
cur_dist = new_dist;
; i <= city_count; ++i)
cur_sol[i] = new_sol[i];
}
else{
; i <= city_count; ++i)
new_sol[i] = cur_sol[i];
}
}
}
T = T*alpha;
}
}
void put_help()
{
puts("Usage: SA data_source <param_option> ");
puts("example: SA C:\\Users\\Administrator\\Desktop\\data.txt");
puts("Format of data file: the coordinates of the cities.");
puts("example:");
puts(" 1 2");
puts(" 3 1");
puts(" 5 6");
puts("Set of param: ");
puts(" -T0 <double> Start temperature");
puts(" -Tf <double> Final temperature");
puts(" -a <double> Cooling coefficient (0 < a < 1)");
puts(" -M <int> Length of Markov chain");
}
void disp_result()
{
printf("最短路径长度为: %lf\n",best_dist);
printf("最优解为: ");
; i <= city_count; ++i){
printf("%d ", best_sol[i]);
}
puts("");
printf("求解用时为: %lf s\n", static_cast<double>(ed - st)/CLOCKS_PER_SEC);
}
main.cpp
用法根据帮助了解
-T0 设置初始温度
-Tf 设置结束温度
-a 设置降温系数
-M 设置Markov链长度
数据文件的路径为绝对路径，内容格式为每行对应记录一个城市的坐标（坐标为二维欧氏空间）
matlab代码：
clear
clc
a = 0.99; % 温度衰减函数的参数
t0 = ; tf = ; t = t0;
Markov_length = ; % Markov链长度
coordinates = [
];
tic
amount = size(coordinates,); % 城市的数目
% 通过向量化的方法计算距离矩阵
dist_matrix = zeros(amount, amount);
coor_x_tmp1 = coordinates(:,) * ones(,amount);
coor_x_tmp2 = coor_x_tmp1';
coor_y_tmp1 = coordinates(:,) * ones(,amount);
coor_y_tmp2 = coor_y_tmp1';
dist_matrix = sqrt((coor_x_tmp1-coor_x_tmp2).^ + ...
(coor_y_tmp1-coor_y_tmp2).^);
sol_new = :amount; % 产生初始解
% sol_new是每次产生的新解；sol_current是当前解；sol_best是冷却中的最好解；
E_current = inf;E_best = inf; % E_current是当前解对应的回路距离；
% E_new是新解的回路距离；
% E_best是最优解的
sol_current = sol_new; sol_best = sol_new;
p = ;
while t>=tf
:Markov_length % Markov链长度
% 产生随机扰动
if (rand < 0.5) % 随机决定是进行两交换还是三交换
% 两交换
ind1 = ; ind2 = ;
while (ind1 == ind2)
ind1 = ceil(rand.*amount);
ind2 = ceil(rand.*amount);
end
tmp1 = sol_new(ind1);
sol_new(ind1) = sol_new(ind2);
sol_new(ind2) = tmp1;
else
% 三交换
ind1 = ; ind2 = ; ind3 = ;
while (ind1 == ind2) || (ind1 == ind3) ...
|| (ind2 == ind3) || (abs(ind1-ind2) == )
ind1 = ceil(rand.*amount);
ind2 = ceil(rand.*amount);
ind3 = ceil(rand.*amount);
end
tmp1 = ind1;tmp2 = ind2;tmp3 = ind3;
% 确保ind1 < ind2 < ind3
if (ind1 < ind2) && (ind2 < ind3)
;
elseif (ind1 < ind3) && (ind3 < ind2)
ind2 = tmp3;ind3 = tmp2;
elseif (ind2 < ind1) && (ind1 < ind3)
ind1 = tmp2;ind2 = tmp1;
elseif (ind2 < ind3) && (ind3 < ind1)
ind1 = tmp2;ind2 = tmp3; ind3 = tmp1;
elseif (ind3 < ind1) && (ind1 < ind2)
ind1 = tmp3;ind2 = tmp1; ind3 = tmp2;
elseif (ind3 < ind2) && (ind2 < ind1)
ind1 = tmp3;ind2 = tmp2; ind3 = tmp1;
end
tmplist1 = sol_new((ind1+):(ind2-));
sol_new((ind1+):(ind1+ind3-ind2+)) = ...
sol_new((ind2):(ind3));
sol_new((ind1+ind3-ind2+):ind3) = ...
tmplist1;
end
%检查是否满足约束
% 计算目标函数值（即内能）
E_new = ;
: (amount-)
E_new = E_new + ...
dist_matrix(sol_new(i),sol_new(i+));
end
% 再算上从最后一个城市到第一个城市的距离
E_new = E_new + ...
dist_matrix(sol_new(amount),sol_new());
if E_new < E_current
E_current = E_new;
sol_current = sol_new;
if E_new < E_best
% 把冷却过程中最好的解保存下来
E_best = E_new;
sol_best = sol_new;
end
else
% 若新解的目标函数值小于当前解的，
% 则仅以一定概率接受新解
if rand < exp(-(E_new-E_current)./t)
E_current = E_new;
sol_current = sol_new;
else
sol_new = sol_current;
end
end
end
t=t.*a; % 控制参数t（温度）减少为原来的a倍
end
toc
disp('最优解为：')
disp(sol_best)
disp('最短距离：')
disp(E_best)
tsp.m
铛铛铛铛，接下来比较两个代码的运行时间。我们以以下31个城市的坐标为测试数据：
1304 2312
3639 1315
4177 2244
3712 1399
3488 1535
3326 1556
3238 1229
4196 1004
4312 790
4386 570
3007 1970
2562 1756
2788 1491
2381 1676
1332 695
3715 1678
3918 2179
4061 2370
3780 2212
3676 2578
4029 2838
4263 2931
3429 1908
3507 2367
3394 2643
3439 3201
2935 3240
3140 3550
2545 2357
2778 2826
2370 2975
先是c语言同学的测试结果
接下来是matlab同学的：
从结果上看，c代码的运行速度远快过matlab代码。
由于matlab代码里考虑了二交换和三交换，在c代码里我偷工减料只写了二交换，所以在最后对对比结果会有影响，但是这个影响微乎其微，我就忽略不计了（。。）
在以上各三次的结果中matlab代码求得的最短距离可能优于c代码的结果，但这并不说明，matlab写的精度较高，毕竟模拟退火是个概率算法。而且用跑一个matlab代码的时间去跑c代码，足以求出更优解吧。
就说到这里了，想到那几天在机房没日没夜跑结果，我去厕所哭会儿。