前几天在做孔群加工问题，各种假设到最后就是求解旅行商问题了，因为原本就有matlab代码模板所以当时就改了城市坐标直接用了，发现运行速度惨不忍睹，最后用上了两个队友的电脑一起跑。这次模拟结束后在想用c语言来实现的话应该可以提高不少效率。关于模拟退火和旅行商问题的介绍我就不赘述了，网上各路大神说的都很详细，我下面就把c语言和matlab代码先附上。
c语言：
1 #ifndef _OPTION_H 2 #define _OPTION_H 3 /* 4 * T0 表示 初始温度 5 * Tf 表示 结束时的温度 6 * alpha 表示 温度衰减系数 7 * T 表示 当前温度 8 * Markov_length 表示 Markov链长度 9 */ 10 double T0 = 100; 11 double alpha = 0.99; 12 double Tf = 3; 13 int Markov_length = 10000; 14 double T = T0; 15 #endif // !_OPTION_H
OPTION.h
1 #ifndef _CITY_H
2 #define _CITY_H
3 struct CITY
4 {
5 int id;
6 double x;
7 double y;
8 };
9 #endif
CITY.h
1 #include <cstdio>
2 #include <cmath>
3 #include <windows.h>
4 #include <ctime>
5 #include <algorithm>
6 #include "OPTION.h"
7 #include "CITY.h"
8 using namespace std;
9 typedef unsigned int UINT;
10 typedef long long LL;
11
12 const double inf = (1 << 30);
13 const double rand_max = 32767;
14 double dist[100][100]; /* 城市之间的距离矩阵 */
15 int new_sol[100]; /* 产生的新解 */
16 int cur_sol[100]; /* 当前解 */
17 int best_sol[100]; /* 冷却中得到的最优解 */
18 double new_dist; /* 新解的路径距离 */
19 double cur_dist; /* 当前解对应路径距离 */
20 double best_dist; /* 冷却过程中得到的最优解 */
21 CITY citys[100]; /* 保存城市坐标和编号信息 */
22 int city_count = 0; /* 保存城市数量 */
23 char *file = NULL; /* 数据文件 */
24 time_t st, ed; /* 记录开始结束时间 */
25
26 void put_help();
27 double init_rand();
28 void input();
29 void cal_dist();
30 void init();
31 void cool();
32 void disp_result();
33
34 int main(int argc, char **argv)
35 {
36 srand(static_cast<int>(time(NULL)));
37 if (argc == 1) {
38 put_help();
39 goto end;
40 }
41 /* get the parameters */
42 file = argv[1];
43 if (argc > 2)
44 {
45 for (int i = 2; i < argc; i += 2) {
46 if (strcmp(argv[i], "-T0") == 0) {
47 T0 = atof(argv[i + 1]);
48 T = T0;
49 }
50 else if (strcmp(argv[i], "-Tf") == 0) {
51 Tf = atof(argv[i + 1]);
52 }
53 else if (strcmp(argv[i], "-a") == 0){
54 alpha = atof(argv[i + 1]);
55 }
56 else if (strcmp(argv[i], "-M") == 0) {
57 Markov_length = atoi(argv[i + 1]);
58 }
59 }
60 }
61 freopen(file, "r", stdin);
62 input();
63 st = clock();
64 cal_dist();
65 init();
66 cool();
67 ed = clock();
68 disp_result();
69 end:
70 return 0;
71 }
72
73 double init_rand()
74 {
75 return rand() / rand_max;
76 }
77
78 void input()
79 {
80 double x, y;
81 while (~scanf("%lf%lf", &x, &y)) {
82 ++city_count;
83 citys[city_count].id = city_count;
84 citys[city_count].x = x;
85 citys[city_count].y = y;
86 }
87 puts("Data input success!");
88 }
89
90 void cal_dist()
91 {
92 for (int i = 1; i <= city_count; ++i) {
93 for (int j = 1; j <= city_count; ++j) {
94 dist[i][j] = sqrt((citys[i].x - citys[j].x)*(citys[i].x - citys[j].x)
95 + (citys[i].y - citys[j].y)*(citys[i].y - citys[j].y));
96 }
97 }
98 }
99
100 void init()
101 {
102 cur_dist = inf;
103 best_dist = inf;
104 for (int i = 1; i <= city_count; ++i) {
105 new_sol[i] = i;
106 cur_sol[i] = i;
107 best_sol[i] = i;
108 }
109 }
110
111 void cool()
112 {
113 while (T >= Tf) {
114 for (int i = 1; i <= Markov_length; ++i) {
115 int ind1 = 0;
116 int ind2 = 0;
117 while (ind1 == ind2 || ind1 == 0 || ind2 == 0 ) {
118 ind1 = static_cast<int>(ceil(city_count*init_rand()));
119 ind2 = static_cast<int>(ceil(city_count*init_rand()));
120 }
121 swap(new_sol[ind1], new_sol[ind2]);
122 /* 检查是否满足约束 */
123 /* 计算距离 */
124 new_dist = 0;
125 for (int i = 1; i <= city_count - 1; ++i) {
126 new_dist = new_dist + dist[new_sol[i]][new_sol[i + 1]];
127 }
128 new_dist += dist[new_sol[city_count]][new_sol[1]];
129 /* 若新解距离小于当前最优距离,则接受新解；否则以一定概率接受新解 */
130 if (new_dist < cur_dist) {
131 cur_dist = new_dist;
132 for (int i = 1; i <= city_count; ++i)
133 cur_sol[i] = new_sol[i];
134 if (new_dist < best_dist) {
135 best_dist = new_dist;
136 for (int i = 1; i <= city_count; ++i)
137 best_sol[i] = new_sol[i];
138 }
139 }
140 else {
141 if (init_rand() < exp(-(new_dist - cur_dist) / T)) {
142 cur_dist = new_dist;
143 for (int i = 1; i <= city_count; ++i)
144 cur_sol[i] = new_sol[i];
145 }
146 else{
147 for (int i = 1; i <= city_count; ++i)
148 new_sol[i] = cur_sol[i];
149 }
150 }
151 }
152 T = T*alpha;
153 }
154 }
155
156 void put_help()
157 {
158 puts("Usage: SA data_source <param_option> ");
159 puts("example: SA C:\\Users\\Administrator\\Desktop\\data.txt");
160 puts("Format of data file: the coordinates of the cities.");
161 puts("example:");
162 puts(" 1 2");
163 puts(" 3 1");
164 puts(" 5 6");
165 puts("Set of param: ");
166 puts(" -T0 <double> Start temperature");
167 puts(" -Tf <double> Final temperature");
168 puts(" -a <double> Cooling coefficient (0 < a < 1)");
169 puts(" -M <int> Length of Markov chain");
170 }
171
172 void disp_result()
173 {
174 printf("最短路径长度为: %lf\n",best_dist);
175 printf("最优解为: ");
176 for (int i = 1; i <= city_count; ++i){
177 printf("%d ", best_sol[i]);
178 }
179 puts("");
180 printf("求解用时为: %lf s\n", static_cast<double>(ed - st)/CLOCKS_PER_SEC);
181 }
main.cpp
用法根据帮助了解
-T0 设置初始温度
-Tf 设置结束温度
-a 设置降温系数
-M 设置Markov链长度
数据文件的路径为绝对路径，内容格式为每行对应记录一个城市的坐标（坐标为二维欧氏空间）
matlab代码：
1 clear 2 clc 3 a = 0.99; % 温度衰减函数的参数 4 t0 = 97; tf = 3; t = t0; 5 Markov_length = 10000; % Markov链长度 6 coordinates = [ 7 1304 2312 8 3639 1315 9 4177 2244 10 3712 1399 11 3488 1535 12 3326 1556 13 3238 1229 14 4196 1004 15 4312 790 16 4386 570 17 3007 1970 18 2562 1756 19 2788 1491 20 2381 1676 21 1332 695 22 3715 1678 23 3918 2179 24 4061 2370 25 3780 2212 26 3676 2578 27 4029 2838 28 4263 2931 29 3429 1908 30 3507 2367 31 3394 2643 32 3439 3201 33 2935 3240 34 3140 3550 35 2545 2357 36 2778 2826 37 2370 2975 38 ]; 39 tic 40 amount = size(coordinates,1); % 城市的数目 41 % 通过向量化的方法计算距离矩阵 42 dist_matrix = zeros(amount, amount); 43 coor_x_tmp1 = coordinates(:,1) * ones(1,amount); 44 coor_x_tmp2 = coor_x_tmp1‘; 45 coor_y_tmp1 = coordinates(:,2) * ones(1,amount); 46 coor_y_tmp2 = coor_y_tmp1‘; 47 dist_matrix = sqrt((coor_x_tmp1-coor_x_tmp2).^2 + ... 48 (coor_y_tmp1-coor_y_tmp2).^2); 49 50 sol_new = 1:amount; % 产生初始解 51 % sol_new是每次产生的新解；sol_current是当前解；sol_best是冷却中的最好解； 52 E_current = inf;E_best = inf; % E_current是当前解对应的回路距离； 53 % E_new是新解的回路距离； 54 % E_best是最优解的 55 sol_current = sol_new; sol_best = sol_new; 56 p = 1; 57 58 while t>=tf 59 for r=1:Markov_length % Markov链长度 60 % 产生随机扰动 61 if (rand < 0.5) % 随机决定是进行两交换还是三交换 62 % 两交换 63 ind1 = 0; ind2 = 0; 64 while (ind1 == ind2) 65 ind1 = ceil(rand.*amount); 66 ind2 = ceil(rand.*amount); 67 end 68 tmp1 = sol_new(ind1); 69 sol_new(ind1) = sol_new(ind2); 70 sol_new(ind2) = tmp1; 71 else 72 % 三交换 73 ind1 = 0; ind2 = 0; ind3 = 0; 74 while (ind1 == ind2) || (ind1 == ind3) ... 75 || (ind2 == ind3) || (abs(ind1-ind2) == 1) 76 ind1 = ceil(rand.*amount); 77 ind2 = ceil(rand.*amount); 78 ind3 = ceil(rand.*amount); 79 end 80 tmp1 = ind1;tmp2 = ind2;tmp3 = ind3; 81 % 确保ind1 < ind2 < ind3 82 if (ind1 < ind2) && (ind2 < ind3) 83 ; 84 elseif (ind1 < ind3) && (ind3 < ind2) 85 ind2 = tmp3;ind3 = tmp2; 86 elseif (ind2 < ind1) && (ind1 < ind3) 87 ind1 = tmp2;ind2 = tmp1; 88 elseif (ind2 < ind3) && (ind3 < ind1) 89 ind1 = tmp2;ind2 = tmp3; ind3 = tmp1; 90 elseif (ind3 < ind1) && (ind1 < ind2) 91 ind1 = tmp3;ind2 = tmp1; ind3 = tmp2; 92 elseif (ind3 < ind2) && (ind2 < ind1) 93 ind1 = tmp3;ind2 = tmp2; ind3 = tmp1; 94 end 95 96 tmplist1 = sol_new((ind1+1):(ind2-1)); 97 sol_new((ind1+1):(ind1+ind3-ind2+1)) = ... 98 sol_new((ind2):(ind3)); 99 sol_new((ind1+ind3-ind2+2):ind3) = ... 100 tmplist1; 101 end 102 103 %检查是否满足约束 104 105 % 计算目标函数值（即内能） 106 E_new = 0; 107 for i = 1 : (amount-1) 108 E_new = E_new + ... 109 dist_matrix(sol_new(i),sol_new(i+1)); 110 end 111 % 再算上从最后一个城市到第一个城市的距离 112 E_new = E_new + ... 113 dist_matrix(sol_new(amount),sol_new(1)); 114 115 if E_new < E_current 116 E_current = E_new; 117 sol_current = sol_new; 118 if E_new < E_best 119 % 把冷却过程中最好的解保存下来 120 E_best = E_new; 121 sol_best = sol_new; 122 end 123 else 124 % 若新解的目标函数值小于当前解的， 125 % 则仅以一定概率接受新解 126 if rand < exp(-(E_new-E_current)./t) 127 E_current = E_new; 128 sol_current = sol_new; 129 else 130 sol_new = sol_current; 131 end 132 end 133 end 134 t=t.*a; % 控制参数t（温度）减少为原来的a倍 135 end 136 toc 137 disp(‘最优解为：‘) 138 disp(sol_best) 139 disp(‘最短距离：‘) 140 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代码，足以求出更优解吧。
就说到这里了，想到那几天在机房没日没夜跑结果，我去厕所哭会儿。
时间： 2024-10-25 03:35:56