什么?我的计算机跑不动了?精确式VS启发式

现在失去的,在未来会以另一种方式归来。
——短短的路走走停停

大家,好·久·不·见·啦~!
在为期两个月的C++速成班结束后,魔鬼的老师布置了大作业,幸运的是这次的大作业可以选择自己找题目。正好小编最近在学一些启发式算法,所以心地善良的小编就想呀,如果能把作业和推文和学习结合在一起,“learning+working+sharing”的形式岂不是很好吗~所以这次给大家带来了大作业的分享——启发式算法和精确式算法的时间消耗比较。
差不多就开始吧!

启发式算法及TSP问题介绍

开始这次小实验之前我们起码要知道什么是启发式算法。
在一些寻找最优解问题中,传统的能够求出最优解的精确式算法(如分支界定、动态规划等)花费的时空过大。启发式算法(heuristic algorithm)是一种能以更快速度找出问题近似最优解的方法。启发式算法通过牺牲一些精度来换取较小的时空消耗,可以看作是解决问题的捷径。
其次,这次小实验问题的载体是老熟人TSP问题。这里再简单介绍一下。
TSP问题(Traveling Salesman Problem,旅行商问题),由威廉哈密顿爵士和英国数学家克克曼T.P.Kirkman于19世纪初提出。问题描述如下:
有若干个城市,任何两个城市之间的距离都是确定的,现要求一旅行商从某城市出发必须经过每一个城市且只在一个城市逗留一次,最后回到出发的城市,问如何事先确定一条最短的线路已保证其旅行的费用最少?
如下图所示:

当然,这次我并不打算用讨厌的美国为例。我们的算例采用柏林52城。
测试数据如下:

算法及代码

在之前关于回溯法的推文中我们曾提到TSP问题,所以这次我们选用的精确式算法依旧是回溯法(我爱ctrl+c)。在启发式算法的选择上,因为小编才刚开始学习,所以把可用的都拿出来比较一下,分别是模拟退火法和迭代局部搜索算法。
可以点击链接学习这三种算法:
【算法进阶】用模拟退火(SA, Simulated Annealing)算法解决旅行商问题
【优化算法】迭代局部搜索算法(Iterated local search)探幽(附C++代码及注释)
(前排怒赞老板的推文,真的是一看就懂,亲测)
代码略微分享一下,因为和老板写的不完全相同(偷偷所一句,老板原先的两段代码都有一丢丢小问题)。可以暂时跳过这一部分,在文末下载代码细看。中间会有些许优化处理并没那么好懂。
代码确实很多,三大段,推荐先跳过。
using namespace std;int const N=14; // 城市数量int Y[N],Ybest[N];double f1,f2;int dis[N][N];int count=0;// 中国31个城市坐标double city_pos[N][2] ={{ 565,575 },{ 25,185 },{ 345,750 },{ 945,685 },{ 845,655 },{ 880,660 },{ 25,230 },{ 525,1000 },{ 580,1175 },{ 650,1130 },{ 1605,620 },{ 1220,580 },{ 1465,200 },{ 1530,5 },/*{ 845,680 },{ 725,370 },{ 145,665 },{ 415,635 },{ 510,875 },{ 560,365 },{ 300,465 },{ 520,585 },{ 480,415 },{ 835,625 },{ 975,580 },{ 1215,245 },{ 1320,315 },{ 1250,400 },{ 660,180 },{ 410,250 },{ 420,555 },{ 575,665 },{ 1150,1160 },{ 700,580 },{ 685,595 },{ 685,610 },{ 770,610 },{ 795,645 },{ 720,635 },{ 760,650 },{ 475,960 },{ 95,260 },{ 875,920 },{ 700,500 },{ 555,815 },{ 830,485 },{ 1170,65 },{ 830,610 },{ 605,625 },{ 595,360 },{ 1340,725 },{ 1740,245 }*/};//函数声明double distance(double *,double *); // 计算两个城市距离void swap(int&,int&); //交换函数double testtime();void setdis();void backtrack(int);double distance(double *p1,double *p2) //计算两点间距离{double dist=0;dist=sqrt(pow(p1[0]-p2[0],2)+pow(p1[1]-p2[1],2));return dist;}void setdis(){for(int i=0;ifor (int j=0;jdis[i][j]=dis[j][i]=distance(city_pos[i],city_pos[j]);}void swap(int&a,int&b) //交换 函数{int temp;temp=a;a=b;b=temp;}void backtrack(int t){if (t==N){ //判断边界。if (f2>f1+dis[Y[t-1]][Y[0]]||f2==0){memcpy(Ybest,Y,N*sizeof(int));f2=f1+dis[Y[t-1]][Y[0]];return;}}else{for (int j=t;j//从t开始,相当于第一轮没有交换 {if(f2>f1 || f2==0){swap(Y[t],Y[j]);f1+=dis[Y[t]][Y[t-1]];backtrack(t+1);f1-=dis[Y[t]][Y[t-1]];swap(Y[t],Y[j]);}}}count++;}double testtime(){time_t start,finish;start = clock(); // 程序运行开始计时setdis();for (int i=0;iY[i]=i;backtrack(0);cout<<"回溯算法。测试点数"<"共迭代" <"次,得到的TSP最优距离为"< "路径为"<<endl; for(int j=0;j// 输出最优路径 {if (j==0)cout<elsecout<<"--->"<}finish = clock(); // 程序结束double duration = ((double)(finish-start))/CLOCKS_PER_SEC; // 计算时间cout<<endl;cout<<"程序运行耗时:"<"秒" <<endl;return duration;}int main(){srand((unsigned)time(NULL)); //初始化随机数种子double time=testtime();cout<<"test: time:"<<endl; int a;cin>>a;return 0;}
//使用模拟退火算法(SA)求解TSP问题 By ZLL.HUST//参考自短短的路走走停停 模拟退火代码 参见微信公众号程序猿声using namespace std;double const T0=50000.0; // 初始温度double const T_min=1e-8;double const q=0.98; // 退火系数double const K=1; //公式中的常数Kint const L=1000; // 每个温度时的迭代次数,即链长int const N=52; // 城市数量// 中国31个城市坐标double city_pos[N][2] ={{ 565,575 },{ 25,185 },{ 345,750 },{ 945,685 },{ 845,655 },{ 880,660 },{ 25,230 },{ 525,1000 },{ 580,1175 },{ 650,1130 },{ 1605,620 },{ 1220,580 },{ 1465,200 },{ 1530,5 },{ 845,680 },{ 725,370 },{ 145,665 },{ 415,635 },{ 510,875 },{ 560,365 },{ 300,465 },{ 520,585 },{ 480,415 },{ 835,625 },{ 975,580 },{ 1215,245 },{ 1320,315 },{ 1250,400 },{ 660,180 },{ 410,250 },{ 420,555 },{ 575,665 },{ 1150,1160 },{ 700,580 },{ 685,595 },{ 685,610 },{ 770,610 },{ 795,645 },{ 720,635 },{ 760,650 },{ 475,960 },{ 95,260 },{ 875,920 },{ 700,500 },{ 555,815 },{ 830,485 },{ 1170,65 },{ 830,610 },{ 605,625 },{ 595,360 },{ 1340,725 },{ 1740,245 }};//函数声明double distance(double *,double *); // 计算两个城市距离double path_len(int *); // 计算路径长度double SA(); //Simulated Annealing模拟退火void swap(int&,int&); //交换函数void create_new(int *); // 产生新解void testtime(double *,double *,int,int *);double distance(double *p1,double *p2) //计算两点间距离{double dis=0;dis=sqrt(pow(p1[0]-p2[0],2)+pow(p1[1]-p2[1],2));return dis;}double path_len(int *city_list) //计算路径长度{double path=0;for (int i=0;i<(N-1);i++)path+=distance(city_pos[city_list[i]],city_pos[city_list[i+1]]);//学长的代码少了下面一行,没有回到原城市,所以写错了哦~!path+=distance(city_pos[city_list[0]],city_pos[city_list[N-1]]) ;return path;}void swap(int&a,int&b) //交换 函数{int temp;temp=a;a=b;b=temp;}void create_new(int *city_list) //产生新解{int i=rand()%N;int j=rand()%N;while(j==i)j=rand()%N;swap(city_list[i],city_list[j]);}double SA(int *bestpath){int count=0;int Ynew[N]; // 用于存放一个解int Ybest[N];for (int i=0;iYbest[i]=i;memcpy(Ynew,Ybest,N*sizeof(int));double f1;double f2=path_len(Ybest);double T=T0;while(T>T_min){for (int i=0;i{create_new(Ynew);f1=path_len(Ynew);double de=f1-f2;if (de<0){memcpy(Ybest,Ynew,N*sizeof(int));f2=f1;}else{double r = ((double)rand())/(RAND_MAX);if(exp(de/(K*T)){memcpy(Ybest,Ynew,N*sizeof(int));f2=f1;}elsememcpy(Ynew,Ybest,N*sizeof(int));}}T*=q;count++;}cout<<"模拟退火算法,初始温度T0="<",降温系数q=" <",每个温度迭代"<"次,共降温"< "次,得到的TSP最优距离为"< "路径为"<<endl; for(int j=0;j// 输出最优路径 {if (j==0)cout<elsecout<<"--->"<}int bestpathlong=path_len(bestpath);if (f2for (int i=0;ibestpath[i]=Ybest[i];return f2;}void testtime(double *time,double *result,int num,int *bestpath){time_t start,finish;start = clock(); // 程序运行开始计时result[num]=SA(bestpath);finish = clock(); // 退火过程结束double duration = ((double)(finish-start))/CLOCKS_PER_SEC; // 计算时间cout<<endl;cout<<"程序运行耗时:"<"秒." <<endl;time[num]=duration;}int main(){srand((unsigned)time(NULL)); //初始化随机数种子double time[10],result[10];int bestpath[N];for(int i=0;i// 初始化总最优路径 bestpath[i]=i;for (int i=0;i<10;i++){cout<<"test"<1<<":"<<endl;testtime(time,result,i,bestpath);}for (int i=0;i<10;i++)cout<<"test"<1<<"\t"<<": time:"<cout<<endl<<"最优路径为:"<<endl;for(int i=0;i// 输出最优路径 {if (i==0)cout<elsecout<<"--->"<}int a;cin>>a;return 0;}
//////////////////////////TSP问题 迭代局部搜索求解代码//代码基于老板代码的框架,并大规模参考//作者:ZLL.HUST//时间:2019-11-28////////////////////////using namespace std;int const CITY_SIZE=52; //城市数量//城市坐标typedef class candidate{public:int x;int y;}city, CITIES;//优化值int Delta[CITY_SIZE][CITY_SIZE];//解决方案typedef class Solution{public:int permutation[CITY_SIZE]; //城市排列int cost; //该排列对应的总路线长度}SOLUTION;// 计算邻域操作优化值int calc_delta(int i, int k, int *tmp, CITIES * cities);//计算两个城市间距离int distance_2city(city c1, city c2);//根据产生的城市序列,计算旅游总距离int cost_total(int * cities_permutation, CITIES * cities);//获取随机城市排列, 用于产生初始解void random_permutation(vector <int> &cities_permutation); //(int * cities_permutation)//颠倒数组中下标begin到end的元素位置, 用于two_opt邻域动作void swap_element(int *p, int begin, int end);//邻域动作 反转index_i <-> index_j 间的元素void two_opt_swap(int *cities_permutation, int *new_cities_permutation, int index_i, int index_j);//本地局部搜索,边界条件 max_no_improvevoid local_search(SOLUTION & best_solution, CITIES * cities, int max_no_improve);//判断接受准则bool AcceptanceCriterion(int *cities_permutation, int *old_cities_permutation, CITIES * cities);//将城市序列分成4块,然后按块重新打乱顺序。//用于扰动函数void double_bridge_move(int *cities_permutation, int * new_cities_permutation,CITIES * cities);//扰动void perturbation(CITIES * cities, SOLUTION &best_solution, SOLUTION ¤t_solution);//迭代搜索void iterated_local_search(SOLUTION & best_solution, CITIES * cities, int max_iterations, int max_no_improve);// 更新Deltavoid Update(int i, int k, int *tmp, CITIES * cities);//Debugvoid test();//berlin52城市坐标,最优解7542好像(这貌似是精确的最优解,我得出的结论是7526)CITIES berlin52[CITY_SIZE] = {{ 565,575 },{ 25,185 },{ 345,750 },{ 945,685 },{ 845,655 },{ 880,660 },{ 25,230 },{ 525,1000 },{ 580,1175 },{ 650,1130 },{ 1605,620 },{ 1220,580 },{ 1465,200 },{ 1530,5 },{ 845,680 },{ 725,370 },{ 145,665 },{ 415,635 },{ 510,875 },{ 560,365 },{ 300,465 },{ 520,585 },{ 480,415 },{ 835,625 },{ 975,580 },{ 1215,245 },{ 1320,315 },{ 1250,400 },{ 660,180 },{ 410,250 },{ 420,555 },{ 575,665 },{ 1150,1160 },{ 700,580 },{ 685,595 },{ 685,610 },{ 770,610 },{ 795,645 },{ 720,635 },{ 760,650 },{ 475,960 },{ 95,260 },{ 875,920 },{ 700,500 },{ 555,815 },{ 830,485 },{ 1170,65 },{ 830,610 },{ 605,625 },{ 595,360 },{ 1340,725 },{ 1740,245 }};int main(){srand((unsigned)time(NULL)); //初始化随机数种子time_t start,finish;start = clock(); // 程序运行开始计时int max_iterations = 600;int max_no_improve = 50;SOLUTION best_solution;iterated_local_search(best_solution, berlin52, max_iterations, max_no_improve);cout << endl<<endl<<"搜索完成! 最优路线总长度 = " << best_solution.cost << endl;cout << "最优访问城市序列如下:" << endl;for (int i = 0; i < CITY_SIZE;i++)cout << best_solution.permutation[i]<<"\t";finish = clock(); // 退火过程结束double duration = ((double)(finish-start))/CLOCKS_PER_SEC; // 计算时间cout<<endl;cout<<"程序运行耗时:"<"秒." <<endl;return 0;}void test(){int jinque[10]={0,4,5,3,9,8,7,2,6,1};cout<endl ;int moni[10]={9,3,5,4,0,1,6,2,7,8};cout<endl ;}//迭代搜索void iterated_local_search(SOLUTION & best_solution, CITIES * cities, int max_iterations, int max_no_improve){SOLUTION current_solution;vector <int> cities_permutation;random_permutation(cities_permutation);for (int i=0;ibest_solution.permutation[i]=cities_permutation[i];best_solution.cost=cost_total(best_solution.permutation, cities);local_search(best_solution, cities, max_no_improve);//cout<<"best_solution.cost"<for(int i=0;i{perturbation(cities, best_solution, current_solution);local_search(current_solution, cities, max_no_improve);if(current_solution.costbest_solution=current_solution;cout << setw(13) << setiosflags(ios::left) <<"迭代搜索 " << i << " 次\t" << "最优解 = " << best_solution.cost << " 当前解 = " << current_solution.cost << endl;}}//计算两个城市间距离int distance_2city(city c1, city c2){double dis=0;dis=sqrt(pow(c1.x-c2.x,2)+pow(c1.y-c2.y,2));return dis;}//根据产生的城市序列,计算旅游总距离int cost_total(int * cities_permutation, CITIES * cities){double path=0;for (int i=0;i<(CITY_SIZE-1);i++)path+=distance_2city(cities[cities_permutation[i]],cities[cities_permutation[i+1]]);path+=distance_2city(cities[cities_permutation[0]],cities[cities_permutation[CITY_SIZE-1]]);return path;}//获取随机城市排列, 用于产生初始解(learning form CSDN)void random_permutation(vector <int> &cities_permutation) //(int * cities_permutation){vector <int> cities_permutation0(CITY_SIZE);for (int i = 0; i < CITY_SIZE; i++)cities_permutation0[i]=i;for (int i=CITY_SIZE;i>0;i--){// 选中的随机下标int index = rand()%i;// 根据选中的下标将原数组选中的元素push到新数组cities_permutation.push_back(cities_permutation0[index]);// 将原数组中选中的元素剔除cities_permutation0.erase(cities_permutation0.begin()+index);}}//颠倒数组中下标begin到end的元素位置void swap_element(int *p, int begin, int end){int temp;while (begin < end){temp = p[begin];p[begin] = p[end];p[end] = temp;begin++;end--;}}//邻域动作 反转index_i <-> index_j 间的元素void two_opt_swap(int *cities_permutation, int *new_cities_permutation, int index_i, int index_j){for (int i = 0; i < CITY_SIZE; i++)new_cities_permutation[i] = cities_permutation[i];swap_element(new_cities_permutation, index_i, index_j);}//本地局部搜索,找出局部最优解,边界条件 max_no_improvevoid local_search(SOLUTION & best_solution, CITIES * cities, int max_no_improve){int count = 0;int inital_cost = best_solution.cost; //初始花费int now_cost = 0;SOLUTION new_solution;for (int i = 0; i < CITY_SIZE - 1; i++)for (int k = i + 1; k < CITY_SIZE; k++) //注意Delta数组中行的数字必定大于列的数字Delta[i][k] = calc_delta(i, k, best_solution.permutation, cities);int deadprotect=0;do{//枚举排列for (int i = 0; i < CITY_SIZE - 1; i++){for (int k = i + 1; k < CITY_SIZE; k++){//邻域动作two_opt_swap(best_solution.permutation, new_solution.permutation, i, k);now_cost = inital_cost + Delta[i][k];new_solution.cost = now_cost;if (new_solution.cost < best_solution.cost){count = 0; //better cost found, so reset.we need to found a best solution in an area(邻域)best_solution = new_solution;inital_cost = best_solution.cost;Update(i, k, best_solution.permutation, cities);}}}count++;deadprotect++;if (deadprotect==4*max_no_improve) break;//cout<} while (count <= max_no_improve);}// 计算邻域操作优化值int calc_delta(int i, int k, int *tmp, CITIES * cities){/*以下计算说明:对于每个方案,翻转以后没必要再次重新计算总距离只需要在翻转的头尾做个小小处理比如:有城市序列 1-2-3-4-5 总距离 = d12 + d23 + d34 + d45 + d51 = A翻转后的序列 1-4-3-2-5 总距离 = d14 + d43 + d32 + d25 + d51 = B由于 dij 与 dji是一样的,所以B也可以表示成 B = A - d12 - d45 + d14 + d25下面的优化就是基于这种原理*///学长这里疯狂分类凑行数吗2333int delta=0;if((i==0)&&(k==CITY_SIZE-1))delta=0;else{int i2=(i-1+CITY_SIZE)%CITY_SIZE;int k2=(k+1)%CITY_SIZE;delta = 0- distance_2city(cities[tmp[i2]], cities[tmp[i]])+ distance_2city(cities[tmp[i2]], cities[tmp[k]])- distance_2city(cities[tmp[k]], cities[tmp[k2]])+ distance_2city(cities[tmp[i]], cities[tmp[k2]]);}return delta;}// 更新Deltavoid Update(int i, int k, int *tmp, CITIES * cities){/*小编注:这是一个计算每次领域动作后路径长度改变值的小优化,在测试中我们暂时忽略这个优化。去重处理,对于Delta数组来说,对于城市序列1-2-3-4-5-6-7-8-9-10,如果对3-5应用了邻域操作2-opt , 事实上对于7-10之间的距离是不需要重复计算的,即后半部分是不变的。 所以用Delta提前预处理一下。当然由于这里的计算本身是O(1) 的,事实上并没有带来时间复杂度的减少(更新操作反而增加了复杂度)如果delta计算 是O(n)的,这种去重操作效果是明显的。if (i && (k != (CITY_SIZE - 1))){ // 如果不是边界,更新(i-1, k + 1)之间的i --; k ++;for (int j = i; j <= k; j ++)for (int l = j + 1; l < CITY_SIZE; l ++)Delta[j][l] = calc_delta(j, l, tmp, cities);for (int j = 0; j < k; j ++){for (int l = i; l <= k; l ++){if (j >= l) continue;Delta[j][l] = calc_delta(j, l, tmp, cities);}}}// else if ((!i)&&(k!=(CITY_SIZE-1))){}// else if (i&&(k==(CITY_SIZE-1))){}//本想在这里加入两个elseif,把两种情况都优化一下,后来想想没必要,其实这都是很小的一部分else*/{for (i = 0; i < CITY_SIZE - 1; i++)for (k = i + 1; k < CITY_SIZE; k++)Delta[i][k] = calc_delta(i, k, tmp, cities);}}//判断接受准则bool AcceptanceCriterion(int *cities_permutation, int *new_cities_permutation, CITIES * cities){int AcceptLimite=500;int c1=cost_total(cities_permutation, cities);int c2=cost_total(new_cities_permutation, cities)-500;if (c1return false;elsereturn true;}//将城市序列分成4块,然后按块重新打乱顺序。//用于扰动函数void double_bridge_move(int *cities_permutation, int * new_cities_permutation,CITIES * cities){int pos1=rand()%(CITY_SIZE/3);int pos2=rand()%(CITY_SIZE/3)+CITY_SIZE/3;int pos3=rand()%(CITY_SIZE/3)+(CITY_SIZE/3)*2;int deadprotect=0;do{int i=0;//写完才发现老板也是反向,哈哈,英雄所见略同啊!for (int j=pos3;j{new_cities_permutation[i]=cities_permutation[j];i++;}for (int j=pos2;j{new_cities_permutation[i]=cities_permutation[j];i++;}for (int j=pos1;j{new_cities_permutation[i]=cities_permutation[j];i++;}for (int j=0;j{new_cities_permutation[i]=cities_permutation[j];i++;}deadprotect++;//cout<//cout<if (deadprotect==5) break;//学长的代码在这里如果没有accept就直接用原数据,所以不需要deadprotect。}while(AcceptanceCriterion(cities_permutation, new_cities_permutation, cities));}//扰动void perturbation(CITIES * cities, SOLUTION &best_solution, SOLUTION ¤t_solution){double_bridge_move(best_solution.permutation, current_solution.permutation,cities);current_solution.cost=cost_total(current_solution.permutation, cities);}
(小错误见注释)
霍霍,万事俱备只欠东风啦!马上开始我们的测试!

测试

先用回溯法跑一跑中国31城的算例:

发现没有输出。为什么呢?

咳咳,并不是程序错了。而是因为回溯法的时间复杂度是指数级的,当n=31时已经是很大的数字了,所以程序一直在跑呀跑呀却出不来。
调整城市的数量,砍到13城左右。

这回总算是跑动啦。
接下来看启发式算法的表现:


(模拟退火,柏林52城问题)
31城、52城都能轻松跑出来。并且耗时都不长,最长不过21秒。
经过小编的多次试验,我们给出三种方法的比较折线图:

(纵坐标单位为秒)

(横坐标表结点数)

从数据中我们可以看出,回溯法的耗时是指数级的,两种启发式算法则是线性增长的。这与理论是吻合的。
从柏林52城的已知解为7542以及迭代局部搜索得出的解7526(应该是由于精度问题导致的差异,我是用整数型表示距离),以及以上的比较结果看来,启发式算法还是比较精确的。
当然我们会发现两种启发式算法所用的时间和精度仍有差别。
时间上的差异是因为人为制定的搜索次数不同,比如模拟退火的退火系数,迭代局部搜索的迭代次数,这些是人为设置的,自然在速度上有差别。比方说我这里设置的模拟退火法明显次数较少,4秒、7秒内就跑完一次搜索(我一次测试了10组);而迭代局部用了7秒、21秒。
至于精度,我认为是算法和代码的差异导致的。在模拟退火中每次都在同一邻域的周围搜索,代码中跳跃的幅度也不大,所以在数据较大时(52城)搜索多次也没有太大的差别,而13城以内比较精确;而迭代局部搜索有一个邻域间跳跃的过程,在此之后进行邻域内搜索最优解,找到整体最优的可能性就大大提高了。
引用老板的一句话结尾:“从上面的过程我们可以看出,模拟退火算法是一种随机算法,它有一定的概率能求得全局最优解,但不一定。用模拟退火算法可以较快速地找出问题的近似最优解。不过,在小编不成熟的眼光看来,人生亦有相似之处。有时候可能放弃眼前短浅的利益,最终才可能获得更好的未来。现在失去的,在未来会以另一种方式归来。”
代码链接: https://pan.baidu.com/s/154PynBJHbfB63eAPl8viaA 提取码: au1f

#
-The End-
文/这里是小编舟寒丶
版/来了这么多新小编不要忘了过去的舟寒丶
码/准备开始猛学一场的舟寒丶
审/帅气老板短短的路走走停停
#
