干货 | 变邻域搜索算法(VNS)求解TSP（附C++详细代码及注释）

  1////////////////////////
  2//TSP问题 变邻域搜索求解代码
  3//基于Berlin52例子求解
  4//作者：infinitor
  5//时间：2018-04-12
  6////////////////////////
  7
  8
  9#include <iostream>
 10#include <cmath>
 11#include <stdlib.h>
 12#include <time.h>
 13#include <vector>
 14#include <windows.h>
 15#include <memory.h>
 16#include <string.h>
 17#include <iomanip>
 18#include <algorithm> 
 19#define DEBUG
 20
 21using namespace std;
 22
 23#define CITY_SIZE 52 //城市数量
 24
 25
 26//城市坐标
 27typedef struct candidate
 28{
 29    int x;
 30    int y;
 31}city, CITIES;
 32
 33//解决方案
 34typedef struct Solution
 35{
 36    int permutation[CITY_SIZE]; //城市排列
 37    int cost;                        //该排列对应的总路线长度
 38}SOLUTION;
 39
 40//城市排列
 41int permutation[CITY_SIZE];
 42//城市坐标数组
 43CITIES cities[CITY_SIZE];
 44
 45
 46//berlin52城市坐标，最优解7542好像
 47CITIES berlin52[CITY_SIZE] =
 48{ 
 49{ 565,575 },{ 25,185 },{ 345,750 },{ 945,685 },{ 845,655 },
 50{ 880,660 },{ 25,230 },{ 525,1000 },{ 580,1175 },{ 650,1130 },{ 1605,620 },
 51{ 1220,580 },{ 1465,200 },{ 1530,5 },{ 845,680 },{ 725,370 },{ 145,665 },
 52{ 415,635 },{ 510,875 },{ 560,365 },{ 300,465 },{ 520,585 },{ 480,415 },
 53{ 835,625 },{ 975,580 },{ 1215,245 },{ 1320,315 },{ 1250,400 },{ 660,180 },
 54{ 410,250 },{ 420,555 },{ 575,665 },{ 1150,1160 },{ 700,580 },{ 685,595 },
 55{ 685,610 },{ 770,610 },{ 795,645 },{ 720,635 },{ 760,650 },{ 475,960 },
 56{ 95,260 },{ 875,920 },{ 700,500 },{ 555,815 },{ 830,485 },{ 1170,65 },
 57{ 830,610 },{ 605,625 },{ 595,360 },{ 1340,725 },{ 1740,245 } 
 58};
 59//优化值
 60int Delta1[CITY_SIZE][CITY_SIZE] = { 0 };
 61
 62
 63//计算两个城市间距离
 64int distance_2city(city c1, city c2)
 65{
 66    int distance = 0;
 67    distance = sqrt((double)((c1.x - c2.x)*(c1.x - c2.x) + (c1.y - c2.y)*(c1.y - c2.y)));
 68
 69    return distance;
 70}
 71
 72//根据产生的城市序列，计算旅游总距离
 73//所谓城市序列，就是城市先后访问的顺序，比如可以先访问ABC，也可以先访问BAC等等
 74//访问顺序不同，那么总路线长度也是不同的
 75//p_perm 城市序列参数
 76int cost_total(int * cities_permutation, CITIES * cities)
 77{
 78    int total_distance = 0;
 79    int c1, c2;
 80    //逛一圈，看看最后的总距离是多少
 81    for (int i = 0; i < CITY_SIZE; i++)
 82    {
 83        c1 = cities_permutation[i];
 84        if (i == CITY_SIZE - 1) //最后一个城市和第一个城市计算距离
 85        {
 86            c2 = cities_permutation[0];
 87        }
 88        else
 89        {
 90            c2 = cities_permutation[i + 1];
 91        }
 92        total_distance += distance_2city(cities[c1], cities[c2]);
 93    }
 94
 95    return total_distance;
 96}
 97
 98//获取随机城市排列
 99void random_permutation(int * cities_permutation)
100{
101    int i, r, temp;
102    for (i = 0; i < CITY_SIZE; i++)
103    {
104        cities_permutation[i] = i; //初始化城市排列，初始按顺序排
105    }
106
107    random_shuffle(cities_permutation, cities_permutation + CITY_SIZE); //随机化排序 
108
109}
110//对应two_opt_swap的去重
111int calc_delta1(int i, int k, int *tmp, CITIES * cities) {
112    int delta = 0;
113    /*
114    以下计算说明：
115    对于每个方案，翻转以后没必要再次重新计算总距离
116    只需要在翻转的头尾做个小小处理
117
118    比如：
119    有城市序列   1-2-3-4-5 总距离 = d12 + d23 + d34 + d45 + d51 = A
120    翻转后的序列 1-4-3-2-5 总距离 = d14 + d43 + d32 + d25 + d51 = B
121    由于 dij 与 dji是一样的，所以B也可以表示成 B = A - d12 - d45 + d14 + d25
122    下面的优化就是基于这种原理
123    */
124    if (i == 0)
125    {
126        if (k == CITY_SIZE - 1)
127        {
128            delta = 0;
129        }
130        else
131        {
132            delta = 0
133                - distance_2city(cities[tmp[k]], cities[tmp[k + 1]])
134                + distance_2city(cities[tmp[i]], cities[tmp[k + 1]])
135                - distance_2city(cities[tmp[CITY_SIZE - 1]], cities[tmp[i]])
136                + distance_2city(cities[tmp[CITY_SIZE - 1]], cities[tmp[k]]);
137        }
138
139    }
140    else
141    {
142        if (k == CITY_SIZE - 1)
143        {
144            delta = 0
145                - distance_2city(cities[tmp[i - 1]], cities[tmp[i]])
146                + distance_2city(cities[tmp[i - 1]], cities[tmp[k]])
147                - distance_2city(cities[tmp[0]], cities[tmp[k]])
148                + distance_2city(cities[tmp[i]], cities[tmp[0]]);
149        }
150        else
151        {
152            delta = 0
153                - distance_2city(cities[tmp[i - 1]], cities[tmp[i]])
154                + distance_2city(cities[tmp[i - 1]], cities[tmp[k]])
155                - distance_2city(cities[tmp[k]], cities[tmp[k + 1]])
156                + distance_2city(cities[tmp[i]], cities[tmp[k + 1]]);
157        }
158    }
159
160    return delta;
161}
162
163
164/*
165去重处理，对于Delta数组来说，对于城市序列1-2-3-4-5-6-7-8-9-10，如果对3-5应用了邻域操作2-opt ， 事实上对于
1667-10之间的翻转是不需要重复计算的。 所以用Delta提前预处理一下。
167
168当然由于这里的计算本身是O（1） 的，事实上并没有带来时间复杂度的减少（更新操作反而增加了复杂度）
169如果delta计算 是O（n）的，这种去重操作效果是明显的。
170*/
171//对应two_opt_swap的去重更新
172void Update1(int i, int k, int *tmp, CITIES * cities, int Delta[CITY_SIZE][CITY_SIZE]) {
173    if (i && k != CITY_SIZE - 1) {
174        i--; k++;
175        for (int j = i; j <= k; j++) {
176            for (int l = j + 1; l < CITY_SIZE; l++) {
177                Delta[j][l] = calc_delta1(j, l, tmp, cities);
178            }
179        }
180
181        for (int j = 0; j < k; j++) {
182            for (int l = i; l <= k; l++) {
183                if (j >= l) continue;
184                Delta[j][l] = calc_delta1(j, l, tmp, cities);
185            }
186        }
187    }// 如果不是边界，更新(i-1, k + 1)之间的 
188    else {
189        for (i = 0; i < CITY_SIZE - 1; i++)
190        {
191            for (k = i + 1; k < CITY_SIZE; k++)
192            {
193                Delta[i][k] = calc_delta1(i, k, tmp, cities);
194            }
195        }
196    }// 边界要特殊更新 
197
198}
199
200
201// two_opt_swap算子 
202void two_opt_swap(int *cities_permutation, int b, int c) 
203{
204    vector<int> v;
205    for (int i = 0; i < b; i++) 
206    {
207        v.push_back(cities_permutation[i]);
208    }
209    for (int i = c; i >= b; i--) 
210    {
211        v.push_back(cities_permutation[i]);
212    }
213    for (int i = c + 1; i < CITY_SIZE; i++) 
214    {
215        v.push_back(cities_permutation[i]);
216    }
217
218    for (int i = 0; i < CITY_SIZE; i++)
219    {
220        cities_permutation[i] = v[i];
221    }
222
223}
224
225//邻域结构1 使用two_opt_swap算子
226void neighborhood_one(SOLUTION & solution, CITIES *cities)
227{
228    int i, k, count = 0;
229    int max_no_improve = 60;
230
231    int inital_cost = solution.cost; //初始花费
232    int now_cost = 0;
233
234    //SOLUTION current_solution = solution;
235
236    for (int i = 0; i < CITY_SIZE - 1; i++)
237    {
238        for (k = i + 1; k < CITY_SIZE; k++)
239        {
240            Delta1[i][k] = calc_delta1(i, k, solution.permutation, cities);
241        }
242    }
243
244    do 
245    {
246        count++;
247        for (i = 0; i < CITY_SIZE - 1; i++)
248        {
249            for (k = i + 1; k < CITY_SIZE; k++)
250            {
251                if (Delta1[i][k] < 0)
252                {
253                    //current_solution = solution;
254                    two_opt_swap(solution.permutation, i, k);
255
256                    now_cost = inital_cost + Delta1[i][k];
257                    solution.cost = now_cost;
258
259                    inital_cost = solution.cost;
260
261                    Update1(i, k, solution.permutation, cities, Delta1);
262
263                    count = 0; //count复位
264
265                }
266
267             }
268          }
269    }while (count <= max_no_improve);
270
271}
272
273//two_h_opt_swap的去重
274int calc_delta2(int i, int k, int *cities_permutation, CITIES * cities)
275{
276    int delta = 0;
277    if (i == 0)
278    {
279        if ( k == i+1)
280        {
281            delta = 0;
282        }
283        else if ( k == CITY_SIZE -1)
284        {
285            delta = 0
286                - distance_2city(cities[cities_permutation[i]], cities[cities_permutation[i + 1]])
287                - distance_2city(cities[cities_permutation[k]], cities[cities_permutation[k - 1]])
288                + distance_2city(cities[cities_permutation[k]], cities[cities_permutation[i+1]])
289                + distance_2city(cities[cities_permutation[k - 1]], cities[cities_permutation[i]]);
290        }
291        else
292        {
293            delta = 0
294                - distance_2city(cities[cities_permutation[i]], cities[cities_permutation[i + 1]])
295                - distance_2city(cities[cities_permutation[k]], cities[cities_permutation[k - 1]])
296                - distance_2city(cities[cities_permutation[k]], cities[cities_permutation[k + 1]])
297                + distance_2city(cities[cities_permutation[k - 1]], cities[cities_permutation[k + 1]])
298                + distance_2city(cities[cities_permutation[i]], cities[cities_permutation[k]])
299                + distance_2city(cities[cities_permutation[k]], cities[cities_permutation[i + 1]]);
300        }
301    }
302    else
303    {
304        if (k == i + 1)
305        {
306            delta = 0;
307        }
308        else if ( k == CITY_SIZE - 1)
309        {
310            delta = 0
311                - distance_2city(cities[cities_permutation[i]], cities[cities_permutation[i + 1]])
312                - distance_2city(cities[cities_permutation[0]], cities[cities_permutation[k]])
313                - distance_2city(cities[cities_permutation[k]], cities[cities_permutation[k-1]])
314                + distance_2city(cities[cities_permutation[k]], cities[cities_permutation[i + 1]])
315                + distance_2city(cities[cities_permutation[k-1]], cities[cities_permutation[0]])
316                + distance_2city(cities[cities_permutation[i]], cities[cities_permutation[k]]);
317        }
318        else
319        {
320            delta = 0
321                - distance_2city(cities[cities_permutation[i]], cities[cities_permutation[i + 1]])
322                - distance_2city(cities[cities_permutation[k]], cities[cities_permutation[k + 1]])
323                - distance_2city(cities[cities_permutation[k]], cities[cities_permutation[k - 1]])
324                + distance_2city(cities[cities_permutation[i]], cities[cities_permutation[k]])
325                + distance_2city(cities[cities_permutation[k]], cities[cities_permutation[i + 1]])
326                + distance_2city(cities[cities_permutation[k - 1]], cities[cities_permutation[k + 1]]);
327
328        }
329    }
330
331    return delta;
332
333}
334
335
336
337//two_h_opt_swap算子
338void two_h_opt_swap(int *cities_permutation, int a, int d) 
339{
340    int n = CITY_SIZE;
341    vector<int> v;
342    v.push_back(cities_permutation[a]);
343    v.push_back(cities_permutation[d]);
344    // i = 1 to account for a already added
345    for (int i = 1; i < n; i++) 
346    {
347        int idx = (a + i) % n;
348        // Ignore d which has been added already
349        if (idx != d) 
350        {
351            v.push_back(cities_permutation[idx]);
352        }
353    }
354
355    for (int i = 0; i < v.size(); i++)
356    {
357        cities_permutation[i] = v[i];
358    }
359
360}
361
362
363//邻域结构2 使用two_h_opt_swap算子
364void neighborhood_two(SOLUTION & solution, CITIES *cities)
365{
366    int i, k, count = 0;
367    int max_no_improve = 60;
368    int inital_cost = solution.cost; //初始花费
369    int now_cost = 0;
370    int delta = 0;
371
372    do
373    {
374        count++;
375        for (i = 0; i < CITY_SIZE - 1; i++)
376        {
377            for (k = i + 1; k < CITY_SIZE; k++)
378            {
379
380                delta = calc_delta2(i, k, solution.permutation, cities);
381
382                if (delta < 0)
383                {
384                    //cout<<"delta = " <<delta<<endl; 
385
386                    two_h_opt_swap(solution.permutation, i, k);
387
388                    now_cost = inital_cost + delta;
389                    solution.cost = now_cost;
390
391                    inital_cost = solution.cost;
392
393                    count = 0; //count复位
394                }
395            }
396        }
397    } while (count <= max_no_improve);
398}
399
400
401//VND
402//best_solution最优解
403//current_solution当前解
404void variable_neighborhood_descent(SOLUTION & solution, CITIES * cities)
405{
406
407    SOLUTION current_solution = solution;
408    int l = 1;
409    cout  <<"=====================VariableNeighborhoodDescent=====================" << endl;
410    while(true)
411    {
412        switch (l)
413        {
414        case 1:
415            neighborhood_one(current_solution, cities);
416            cout << setw(45) << setiosflags(ios::left)  <<"Now in neighborhood_one , current_solution = " << current_solution.cost << setw(10) << setiosflags(ios::left) << "  solution = " << solution.cost << endl;
417            if (current_solution.cost < solution.cost)
418            {
419                solution = current_solution;
420                l = 0;
421            }
422            break;
423        case 2:
424            neighborhood_two(current_solution, cities);
425            cout << setw(45) << setiosflags(ios::left) << "Now in neighborhood_two , current_solution = " << current_solution.cost << setw(10) << setiosflags(ios::left) << "  solution = " << solution.cost << endl;
426            if (current_solution.cost < solution.cost)
427            {
428                solution = current_solution;
429                l = 0;
430            }
431            break;
432
433        default:
434            return;
435        }
436        l++;
437
438    }
439
440}
441
442//将城市序列分成4块，然后按块重新打乱顺序。
443//用于扰动函数
444void double_bridge_move(int * cities_permutation)
445{
446    int pos1 = 1 + rand() % (CITY_SIZE / 4);
447    int pos2 = pos1 + 1 + rand() % (CITY_SIZE / 4);
448    int pos3 = pos2 + 1 + rand() % (CITY_SIZE / 4);
449
450    int i;
451    vector<int> v;
452    //第一块
453    for (i = 0; i < pos1; i++)
454    {
455        v.push_back(cities_permutation[i]);
456    }
457
458    //第二块
459    for (i = pos3; i < CITY_SIZE; i++)
460    {
461        v.push_back(cities_permutation[i]);
462    }
463    //第三块
464    for (i = pos2; i < pos3; i++)
465    {
466        v.push_back(cities_permutation[i]);
467    }
468
469    //第四块
470    for (i = pos1; i < pos2; i++)
471    {
472        v.push_back(cities_permutation[i]);
473    }
474
475
476    for (i = 0; i < (int)v.size(); i++)
477    {
478        cities_permutation[i] = v[i];
479    }
480
481
482}
483
484//抖动
485void shaking(SOLUTION &solution, CITIES *cities)
486{
487    double_bridge_move(solution.permutation);
488    solution.cost = cost_total(solution.permutation, cities);
489}
490
491
492void variable_neighborhood_search(SOLUTION & best_solution, CITIES * cities)
493{
494
495    int max_iterations = 5;
496
497    int count = 0, it = 0;
498
499    SOLUTION current_solution = best_solution;
500
501    //算法开始
502    do 
503    {
504        cout << endl << "\t\tAlgorithm VNS iterated  " << it+1 << "  times" << endl;
505        count++;
506        it++;
507        shaking(current_solution, cities);
508
509        variable_neighborhood_descent(current_solution, cities); 
510
511        if (current_solution.cost < best_solution.cost)
512        {
513            best_solution = current_solution;
514            count = 0;
515        }
516
517        cout << "\t\t全局best_solution = " << best_solution.cost << endl;
518
519    } while (count <= max_iterations);
520
521
522}
523
524
525int main()
526{
527
528    srand((unsigned) time(0));
529
530    SOLUTION best_solution;
531
532    random_permutation(best_solution.permutation);
533    best_solution.cost = cost_total(best_solution.permutation, berlin52);
534
535    cout << "初始总路线长度 = " << best_solution.cost << endl;
536
537    variable_neighborhood_search(best_solution, berlin52);
538
539    cout << endl << endl << "搜索完成！ 最优路线总长度 = " << best_solution.cost << endl;
540    cout << "最优访问城市序列如下：" << endl;
541    for (int i = 0; i < CITY_SIZE; i++)
542    {
543        cout << setw(4) << setiosflags(ios::left) << best_solution.permutation[i];
544    }
545
546    cout << endl << endl;
547
548    return 0;
549}