模拟退火算法求解旅行商问题（附c和matlab源代码）

时间：2015-08-19 00:10:08 阅读：251 评论：0 收藏：0 [点我收藏+]

标签：

前几天在做孔群加工问题，各种假设到最后就是求解旅行商问题了，因为原本就有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