模擬退火算法研究

本文轉載自查看原文 2018-01-19 20:53 1020 Research

模擬退火算法的物理背景：

模擬退火算法的核心思想與熱力學的原理極為相似，在高溫下，液體的大量分子彼此之間進行相對移動，如果液體慢慢冷卻，原子的可動性就會消失，原子自行排列成行，形成一個純凈的晶體．如果溫度冷卻迅速，那它不一定能達到這個狀態，這一個過程的本質在於，慢慢冷卻，使原子在喪失可動性之前重新排列，達到低能狀態的必要條件．

模擬退火算法的描述：

模擬退火算法實質是一種貪心算法，但和貪心算法不同的是，模擬退火算法在搜索過程中加入了隨機的因素，它在搜索的過程中，會以一定的概率來接受比當前解要更差的解，因此有可能會跳出這個局部最優解，找到全局的最優值．

如圖，假設從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 }

驗證數據：

執行過程

答案大概在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左右

總結：

模擬退火算法計算過程簡單，引入概率跳出局部最優解，全局性較強．但在執行過程中，算法的執行時間較長，相關參數的變化對算法性能影響大，例如，如果在某一溫度下充分搜索，很容易陷入到局部最優解中，並且收斂速度慢，但如果不充分搜索，那么很難的到全局最優解．所以，約束條件的選擇，參數控制都非常關鍵，還有待進一步研究．

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 模擬退火算法模擬退火算法詳解模擬退火算法和遺傳算法模擬退火算法(1) - Python實現模擬退火算法學習筆記模擬退火算法及其Matlab實現 Matlab 模擬退火算法模型代碼模擬退火模擬退火 [轉]爬山算法和模擬退火算法簡介