模擬退火算法研究
阿新 • • 發佈：2017-12-20
main bre 之前 sca 狀態 ios 最優 流程圖 概率
模擬退火算法的物理背景：
模擬退火算法的核心思想與熱力學的原理極為相似，在高溫下，液體的大量分子彼此之間進行相對移動，如果液體慢慢冷卻，原子的可動性就會消失，原子自行排列成行，形成一個純凈的晶體． 如果溫度冷卻迅速，那它不一定能達到這個狀態，這一個過程的本質在於，慢慢冷卻，使原子在喪失可動性之前重新排列，達到低能狀態的必要條件．
模擬退火算法的描述：
模擬退火算法實質是一種貪心算法，但和貪心算法不同的是，模擬退火算法在搜索過程中加入了隨機的因素，它在搜索的過程中，會以一定的概率來接受比當前解要更差的解，因此有可能會跳出這個局部最優解，找到全局的最優值．
如圖，假設從A點出發，搜索過程中找到了局部最優解B,模擬退火算法會以一定的概率接受比B更差的解，從而找到全局最優解C
模擬退火算法流程圖
模擬退火算法主要分三個步驟：
１．在當前解的基礎上產生一個新 解
２．利用Metropolis準則(即以e(- ?E/T)的概率接受最優解)，判 斷新解是否被接受．
３．當新解確定被接受時，用新解 代替當前解，在此基礎上進行 下一輪實驗，當新解被判定舍 棄時，在原解的基礎上進行下 一輪實驗
模擬退火算法偽代碼：
模擬退火算法解決實際問題：
１．TSP問題：
旅行商問題(Traveling Salesman Problem,TSP),是數學領域一個著名問題之一，假設有一個旅行商人要拜訪n個城市，他必須選擇所要走的路徑，路徑的限制是每個城市只能拜訪一次，而且最後要回到原來出發的城市，路徑的選擇目標是要求的路徑的路程為所有路徑中的最小值 .
執行代碼
1 #include <stdio.h> 2 #include <stdlib.h> 3 #include <math.h> 4 #include <time.h> 5 #include <algorithm> 6 #include <iostream> 7 using namespace std; 8 9 int N; //城市的數量 10 const double T=3000; //初始溫度 11 const double k=0.98; //溫度下降比例12 const double mint=1e-8; //溫度的極值 13 14 const int Oloop=1000; //溫度控制的最小值為t*0.98^20 15 const int Iloop=15000; //尋找一定溫度下的最大值 16 const int lim=10000; //概率選擇次數的最大值 17 18 double dist[100][100]; //任意兩座城市之間的距離 19 20 struct Path 21 { 22 int city[100]; 23 double len; //路徑的總長度 24 }; 25 26 27 struct Point //城市點的坐標 28 { 29 double x; 30 double y; 31 }; 32 33 Path bestpath; //記錄最優路徑 34 Path currpath; //記錄當前的路徑 35 Point p[100]; //各點的坐標 36 37 void input() 38 { 39 scanf("%d",&N); 40 for(int i=0;i<N;i++) 41 { 42 scanf("%lf%lf",&p[i].x,&p[i].y); 43 } 44 } 45 46 //計算任意兩點之間的距離 47 void compu() 48 { 49 for(int i=0;i<N;i++) 50 { 51 for(int j=i+1;j<N;j++) 52 { 53 dist[i][j]=dist[j][i]=sqrt((p[i].x-p[j].x)*(p[i].x-p[j].x)+(p[i].y-p[j].y)*(p[i].y-p[j].y)); 54 } 55 } 56 } 57 58 //初始化最優路徑 59 void init() 60 { 61 bestpath.len=0; 62 bestpath.city[0]=0; 63 for(int i=1;i<N;i++) 64 { 65 bestpath.city[i]=i; 66 bestpath.len=bestpath.len+dist[i-1][i]; 67 } 68 bestpath.len=bestpath.len+dist[0][N-1]; 69 } 70 71 //隨機選擇目的優化路徑，此處有兩種方式，一是隨機選取兩個節點交換，二是隨機排列 72 Path choose(Path newpath) 73 { 74 int city1,city2; 75 city1=(int)rand()%N; 76 city2=(int)rand()%N; 77 while(city1==city2) 78 { 79 city1=(int)rand()%N; 80 city2=(int)rand()%N; 81 } 82 swap(newpath.city[city1],newpath.city[city2]); 83 newpath.len=0; 84 for(int i=1;i<N;i++) 85 { 86 newpath.len=newpath.len+dist[newpath.city[i-1]][newpath.city[i]]; 87 } 88 newpath.len=newpath.len+dist[newpath.city[0]][newpath.city[N-1]]; 89 return newpath; 90 } 91 92 //模擬退火算法核心代碼 93 Path sa() 94 { 95 srand((unsigned)(time(NULL)));//時刻更新隨機模板 96 Path newpath=bestpath; 97 Path currpath=bestpath; 98 int pl=0; //以一定概率接受此路徑的次數 99 int pf=0; 100 double t2=T; 101 double t; 102 while(true){ 103 for(int i=0;i<Iloop;i++) 104 { 105 newpath=choose(currpath); 106 t=newpath.len-currpath.len; 107 if(t<0) 108 { 109 currpath=newpath; 110 pl=0; 111 pf=0; 112 113 } 114 else 115 { 116 double ran=rand()/(RAND_MAX+1.0); 117 double exx=-t*1.0/t2*1.0; 118 double ex=exp(exx); 119 if(ran<ex) 120 { 121 currpath=newpath; 122 } 123 pl++; 124 } 125 if(pl>lim){ 126 pf++; 127 break; 128 } 129 } 130 if(currpath.len<bestpath.len) 131 { 132 bestpath=currpath; 133 } 134 if(pf>Oloop||t2<mint) 135 { 136 break; 137 } 138 t2=t2*k; 139 } 140 } 141 142 void Print() 143 { 144 for(int i=0;i<N;i++) 145 { 146 printf("%d ",bestpath.city[i]); 147 } 148 printf("length：%lf\n",bestpath.len); 149 } 150 151 int main() 152 { 153 for(int i=0;i<50;i++){ 154 freopen("in.txt","r",stdin); 155 input(); 156 compu(); 157 init(); 158 sa(); 159 Print(); 160 } 161 return 0; 162 }
驗證數據：
1 27 2 41 94 3 37 84 4 53 67 5 25 62 6 7 64 7 2 99 8 68 58 9 71 44 10 54 62 11 83 69 12 64 60 13 18 54 14 22 60 15 83 46 16 91 38 17 25 38 18 24 42 19 58 69 20 71 71 21 74 78 22 87 76 23 18 40 24 13 40 25 82 7 26 62 32 27 58 35 28 45 21
執行過程
答案大概在401左右
２．求函數f(x,y)=5cos(xy)+xy+y^3的最小值，其中x的取值範圍是(-5,5),y的取值範圍是(-5,5)
執行代碼：
1 #include <stdio.h> 2 #include <math.h> 3 #include <time.h> 4 #include <stdlib.h> 5 #include <algorithm> 6 #include <iostream> 7 8 using namespace std; 9 10 const double MAXx=5.0; //x的最大值 11 const double MINx=-5.0; //x的最小值 12 const double MAXy=5.0; //y的最大值 13 const double MINy=-5.0; //y的最小值 14 const double T=100; //初始溫度 15 const double S=0.5; //每次的變化量 16 const double K=0.99; //溫度遞減率 17 const double ex=1e-9; 18 double Bestx; 19 double Besty; 20 21 double functionf(double x1,double y1) //所要求解的函數 22 { 23 return 5*cos(x1*y1)+x1*y1+y1*y1*y1; 24 } 25 26 void sa() 27 { 28 srand((unsigned)(time(NULL))); 29 double currx,curry; 30 double newx,newy; 31 double Bestfun,currfun,newfun; 32 double t=T; 33 double de; 34 Bestx=MINx+(MAXx-MINx)*rand()/RAND_MAX; 35 Besty=MINy+(MAXx-MINy)*rand()/RAND_MAX; 36 Bestfun=functionf(Bestx,Besty); 37 currx=Bestx; 38 curry=Besty; 39 while(t>0.00001) 40 { int P=0; 41 for(int i=0;i<10;i++) 42 { 43 int p=0; 44 while(p==0) 45 { 46 newx=currx+S*(MINx+(MAXx-MINx)*rand()/RAND_MAX); 47 if(newx>=MINx&&newx<=MAXx) 48 { 49 p=1; 50 } 51 } 52 de=functionf(currx,curry)-functionf(newx,newy); 53 double dexp=exp(de/t); 54 if(de>0) 55 { 56 currx=newx; 57 //P++; 58 } 59 else 60 { 61 if(dexp>(rand()/RAND_MAX+1)) 62 { 63 currx=newx; 64 //P++; 65 } 66 } 67 for(int j=0;j<10;j++) 68 { 69 p=0; 70 while(p==0){ 71 newy=curry+S*(MINy+(MAXy-MINy)*rand()/RAND_MAX); 72 if(newy>=MINy&&newy<=MAXy) 73 { 74 p=1; 75 } 76 } 77 de=functionf(currx,curry)-functionf(newx,newy); 78 double dexp=exp(de/t); 79 if(de>0) 80 { 81 curry=newy; 82 } 83 else 84 { 85 if(dexp>rand()/RAND_MAX+1) 86 { 87 curry=newy; 88 } 89 } 90 } 91 92 93 94 // double y2=curry; 95 // for(int j=0;j<10;j++) 96 // { 97 // p=0; 98 // while(p==0){ 99 // newy=y2+S*(MINy+(MAXy-MINy)*rand()/RAND_MAX); 100 // if(newy>=MINy&&newy<=MAXy) 101 // { 102 // p=1; 103 // } 104 // } 105 // de=functionf(currx,y2)-functionf(newx,newy); 106 // double dexp=exp(de/t); 107 // if(de>0) 108 // { 109 // y2=newy; 110 // } 111 // else 112 // { 113 // if(dexp>rand()/RAND_MAX+1) 114 // { 115 // y2=newy; 116 // } 117 // } 118 // } 119 // 120 // if(functionf(currx,y2)<functionf(currx,curry)) 121 // { 122 // curry=y2; 123 // } 124 // if(P>300){ 125 // P=0; 126 // break; 127 // } 128 129 130 } 131 if(functionf(currx,curry)<functionf(Bestx,Besty)) 132 { 133 Bestx=currx; 134 Besty=curry; 135 } 136 t=t*K; 137 } 138 } 139 140 int main() 141 { 142 for(int i=1;i<205;i++) 143 { 144 sa(); 145 printf("%f ",functionf(Bestx,Besty)); 146 if(i%6==0) 147 printf("\n"); 148 } 149 printf("\n"); 150 }
執行過程
答案在-151.1左右
總結：
模擬退火算法計算過程簡單，引入概率跳出局部最優解，全局性較強． 但在執行過程中，算法的執行時間較長，相關參數的變化對算法性能影響大，例如，如果在某一溫度下充分搜索，很容易陷入到局部最優解中，並且收斂速度慢，但如果不充分搜索，那麽很難的到全局最優解． 所以，約束條件的選擇，參數控制都非常關鍵，有待進一步研究．
模擬退火算法研究