爬山演算法&模擬退火
優化演算法入門系列文章目錄（更新中）：
2.遺傳演算法
一. 爬山演算法 ( Hill Climbing )
介紹模擬退火前，先介紹爬山演算法。爬山演算法是一種簡單的貪心搜尋演算法，該演算法每次從當前解的臨近解空間中選擇一個最優解作為當前解，直到達到一個區域性最優解。
爬山演算法實現很簡單，其主要缺點是會陷入區域性最優解，而不一定能搜尋到全域性最優解。如圖1所示：假設C點為當前解，爬山演算法搜尋到A點這個區域性最優解就會停止搜尋，因為在A點無論向那個方向小幅度移動都不能得到更優的解。
圖1
version1：
二. 模擬退火(SA,Simulated Annealing)思想
爬山法是完完全全的貪心法，每次都鼠目寸光的選擇一個當前最優解，因此只能搜尋到區域性的最優值。模擬退火其實也是一種貪心演算法，但是它的搜尋過程引入了隨機因素。模擬退火演算法以一定的概率來接受一個比當前解要差的解，因此有可能會跳出這個區域性的最優解，達到全域性的最優解。以圖1為例，模擬退火演算法在搜尋到區域性最優解A後，會以一定的概率接受到E的移動。也許經過幾次這樣的不是區域性最優的移動後會到達D點，於是就跳出了局部最大值A。
模擬退火演算法描述：
若J( Y(i+1) )>= J( Y(i) ) (即移動後得到更優解)，則總是接受該移動
若J( Y(i+1) )< J( Y(i) ) (即移動後的解比當前解要差)，則以一定的概率接受移動，而且這個概率隨著時間推移逐漸降低（逐漸降低才能趨向穩定）
這裡的“一定的概率”的計算參考了金屬冶煉的退火過程，這也是模擬退火演算法名稱的由來。
根據熱力學的原理，在溫度為T時，出現能量差為dE的降溫的概率為P(dE)，表示為：
P(dE) = exp( dE/(kT) )
其中k是一個常數，exp表示自然指數，且dE<0。這條公式說白了就是：溫度越高，出現一次能量差為dE的降溫的概率就越大；溫度越低，則出現降溫的概率就越小。又由於dE總是小於0（否則就不叫退火了），因此dE/kT < 0 ，所以P(dE)的函式取值範圍是(0,1) 。
隨著溫度T的降低，P(dE)會逐漸降低。
我們將一次向較差解的移動看做一次溫度跳變過程，我們以概率P(dE)來接受這樣的移動。
關於爬山演算法與模擬退火，有一個有趣的比喻：
爬山演算法：兔子朝著比現在高的地方跳去。它找到了不遠處的最高山峰。但是這座山不一定是珠穆朗瑪峰。這就是爬山演算法，它不能保證區域性最優值就是全域性最優值。
模擬退火：兔子喝醉了。它隨機地跳了很長時間。這期間，它可能走向高處，也可能踏入平地。但是，它漸漸清醒了並朝最高方向跳去。這就是模擬退火。
下面給出模擬退火的偽程式碼表示。
三. 模擬退火演算法虛擬碼
/*
* J(y)：在狀態y時的評價函式值
* Y(i)：表示當前狀態
* Y(i+1)：表示新的狀態
* r： 用於控制降溫的快慢
* T： 系統的溫度，系統初始應該要處於一個高溫的狀態
* T_min ：溫度的下限，若溫度T達到T_min，則停止搜尋
*/
while( T > T_min )
{
dE = J( Y(i+1) ) - J( Y(i) ) ;
if ( dE >=0 ) //表達移動後得到更優解，則總是接受移動
Y(i+1) = Y(i) ; //接受從Y(i)到Y(i+1)的移動
else
{
// 函式exp( dE/T )的取值範圍是(0,1) ，dE/T越大，則exp( dE/T )也
if ( exp( dE/T ) > random( 0 , 1 ) )
Y(i+1) = Y(i) ; //接受從Y(i)到Y(i+1)的移動
}
T = r * T ; //降溫退火 ，0<r<1 。r越大，降溫越慢；r越小，降溫越快
/*
* 若r過大，則搜尋到全域性最優解的可能會較高，但搜尋的過程也就較長。若r過小，則搜尋的過程會很快，但最終可能會達到一個區域性最優值
*/
i ++ ;
}version2:
模擬退火演算法是用來求解最優化問題的演算法。比如著名的TSP問題，函式最大值最小值問題等等。接下來將以如下幾個方面來詳細介紹模擬退火演算法。
1. 模擬退火演算法認識
爬山演算法也是一個用來求解最優化問題的演算法，每次都向著當前上升最快的方向往上爬，但是初始化不同可能
會得到不同的區域性最優值，模擬退火演算法就可能跳出這種區域性最優解的限制。模擬退火演算法是模擬熱力學系統
中的退火過程。在退火過程中是將目標函式作為能量函式。大致過程如下
初始高溫 => 溫度緩慢下降=> 終止在低溫 （這時能量函式達到最小，目標函式最小）
在熱力學中的退火過程大致是變溫物體緩慢降溫而達到分子之間能量最低的狀態。設熱力學系統S中有有限個且
離散的n個狀態，狀態
模擬退火演算法也是貪心演算法，但是在其過程中引入了隨機因素，以一定的概率接受一個比當前解要差的解，並且
這個概率隨著時間的推移而逐漸降低。
2. 模擬退火演算法描述
若
若
逐漸降低。這個概率表示為
由於是退火過程，所以dE < 0，這個公式說明了溫度越高出現一次能量差為dE的降溫概率就越大，溫度越低，
出現降溫的概率就越小，由於dE總是小於0，所以P(dE)取值在0到1之間。偽碼如下
模擬退火樣例1：
費馬點問題求解
題意：給n個點，找出一個點，使這個點到其他所有點的距離之和最小，也就是求費馬點。
由於現在還未能掌握模擬退火的具體方法，不知怎麼取隨機數，看了好多個版本的程式碼還是未能想明白，所以先貼程式碼先，以後再思考。
程式碼1：
這個方法不太具有普遍性，有些題目可能不太適用，因為它是一個點朝四個方向按一定方向搜尋，然後再逐步縮小步長。程式碼寫的簡單多了，但這種做法有時候會錯。。
#include <iostream>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#define N 1005
#define eps 1e-8 //搜尋停止條件閥值
#define INF 1e99
#define delta 0.98 //溫度下降速度
#define T 100 //初始溫度
using namespace std;
int dx[4] = {0, 0, -1, 1};
int dy[4] = {-1, 1, 0, 0}; //上下左右四個方向
struct Point
{
double x, y;
};
Point p[N];
double dist(Point A, Point B)
{
return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}
double GetSum(Point p[], int n, Point t)
{
double ans = 0;
while(n--)
ans += dist(p[n], t);
return ans;
}
//其實我覺得這玩意兒根本不叫模擬退火
double Search(Point p[], int n)
{
Point s = p[0]; //隨機初始化一個點開始搜尋
double t = T; //初始化溫度
double ans = INF; //初始答案值
while(t > eps)
{
bool flag = 1;
while(flag)
{
flag = 0;
for(int i = 0; i < 4; i++) //上下左右四個方向
{
Point z;
z.x = s.x + dx[i] * t;
z.y = s.y + dy[i] * t;
double tp = GetSum(p, n, z);
if(ans > tp)
{
ans = tp;
s = z;
flag = 1;
}
}
}
t *= delta;
}
return ans;
}
int main()
{
int n;
while(scanf("%d", &n) != EOF)
{
for(int i = 0; i < n; i++)
scanf("%lf %lf", &p[i].x, &p[i].y);
printf("%.0lf\n", Search(p, n));
}
return 0;
}
程式碼2：
這個程式碼符合上面模擬退火的虛擬碼，但是下一個狀態到底應該怎麼取？隨便？我看有些方法是乘以步長，有些乘以角度等等。
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<ctime>
using namespace std;
int n;
double xx,yy,ans,t;
struct point{double x,y;}p[105];
double sqr(double x){return x*x;}
double dis(double x,double y,point p)
{return sqrt(sqr(x-p.x)+sqr(y-p.y));}
double getsum(double x,double y)
{
double tmp=0;
for(int i=1;i<=n;i++)
tmp+=dis(x,y,p[i]);
return tmp;
}
int main()
{
srand(time(0));
while(scanf("%d",&n)!=EOF)
{
xx=yy=0;ans=1e20;t=100000;
for(int i=1;i<=n;i++)
{
scanf("%lf%lf",&p[i].x,&p[i].y);
xx+=p[i].x;yy+=p[i].y;
}
xx/=n;yy/=n;
ans=getsum(xx,yy);
double tmp,x,y;
while(t>0.02)
{
x=y=0;
for(int i=1;i<=n;i++)
{
x+=(p[i].x-xx)/dis(xx,yy,p[i]);
y+=(p[i].y-yy)/dis(xx,yy,p[i]);
}
tmp=getsum(xx+x*t,yy+y*t);
if(tmp<ans)
{ans=tmp;xx+=x*t,yy+=y*t;}
else if(log((tmp-ans)/t)<(rand()%10000)/10000.0)
{ans=tmp;xx+=x*t,yy+=y*t;}
t*=0.9;
}
printf("%.0lf\n",ans);
}
return 0;
}程式碼3：
下面的程式碼也符合虛擬碼的步驟，但是沒有以一定概率接受新的狀態那個程式碼，到底什麼時候需要？什麼時候不需要？
大致步驟如下：
1、隨機選取一個合適的控制條件T作為開始
2、隨機選取P個起始點，作為可行解
3、分別依據內能更新這P個可行解
4、減小控制條件T，直到終止條件
#include<iostream>
#include<cmath>
#include<ctime>
#include<cstdio>
using namespace std;
const int MAXN=100+10;
const int KMEAN=30;
const double INF=1<<30;
const double eps=1e-8;
const double PI=3.141592653;
double x=10000.0,y=10000.0;
int n;
struct Point
{
double x,y;
}myPoint[MAXN],myK[KMEAN];
double dis[KMEAN];
double GetMinDis(double tx,double ty)
{
double temp=0;
for(int i=0;i<n;i++)
{
temp+=sqrt((tx-myPoint[i].x)*(tx-myPoint[i].x)+(ty-myPoint[i].y)*(ty-myPoint[i].y));
}
return temp;
}
int main()
{
int i,j;
double temp;
while(scanf("%d",&n)!=EOF)
{
srand((unsigned)time(NULL));
for(i=0;i<n;i++)
scanf("%lf%lf",&myPoint[i].x,&myPoint[i].y);
for(i=0;i<KMEAN;i++)
{
dis[i]=INF;
myK[i].x=rand()%10001/10000.0*x;
myK[i].y=rand()%10001/10000.0*y;
}
for(i=0;i<KMEAN;i++)
dis[i]=GetMinDis(myK[i].x,myK[i].y);
double theta,delta=10000.0/sqrt(1.0*n);
while(delta>eps)
{
for(i=0;i<KMEAN;i++)
{
double nx=myK[i].x,ny=myK[i].y;
for(j=0;j<KMEAN;j++)
{
double tx,ty;
theta=double(rand()%10001)/10000.0*2*PI;
tx=nx+delta*cos(theta);//可改變方向
ty=ny+delta*sin(theta);
temp=GetMinDis(tx,ty);
if(temp<dis[i])
{
myK[i].x=tx;
myK[i].y=ty;
dis[i]=temp;
}
}
}
delta=0.98*delta;
}
temp=INF;
for(i=0;i<KMEAN;i++)
{
if(dis[i]<temp)
{
temp=dis[i];
}
}
printf("%.0lf\n",temp);
}
return 0;
}這個版本的和程式碼3類似，都是隨機K個可行解，不斷更新可行解，最後再去可行解中的最小值。
//求n邊形的費馬點
//即找到一個點使得這個點到n個點的距離之和最小
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>
const double inf = 1e10;
const double pi = acos(-1.0);
const int Rp = 4;//初始時隨機選擇一些點，不用太多
const int shift = 60;//但是方向一定要多
struct point {
double x,y;
void goto_rand_dir(double key)
{
double d=2*pi*(double)rand()/RAND_MAX;
x+=key*sin(d);
y+=key*cos(d);
}
void Get_Rand_Point(int a,int b)
{
x=rand()%a+1;
y=rand()%b+1;
}
}p[1010],randp[Rp];
double Dis(point a,point b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
double dis[Rp];
int main()
{
int i,j,x,y,k,m;
srand(time(NULL));
while(scanf("%d",&m)!=EOF)
{
for(i=0;i<m;i++)
{
scanf("%lf%lf",&p[i].x,&p[i].y);
}
x=10000;y=10000;
double tmp;
for(i=0;i<Rp;i++)
{
dis[i]=inf;
randp[i].Get_Rand_Point(x,y);
tmp=0;
for(j=0;j<m;j++)
{
tmp+=Dis(randp[i],p[j]);
}
if(tmp<dis[i])
dis[i]=tmp;
}
double key=sqrt(1.0*(x*x+y*y));//初始的步長，要保證初始點從該點出發肯定能包括整個區域
while(key>=0.01)
{
for(i=0;i<Rp;i++)
{
for(j=0;j<shift;j++)
{
point cc=randp[i];
cc.goto_rand_dir(key);
if(cc.x<0||cc.y<0||cc.x>x||cc.y>y) continue;
tmp=0;
for(k=0;k<m;k++)
{
tmp+=Dis(cc,p[k]);
}
if(tmp<dis[i]) //如果從i點出發隨機移動的點比原來的點更優，則接受該移動
{
dis[i]=tmp;
randp[i]=cc;
}
}
}
key=key*0.6;//可靈活調整
}
for(i=k=0;i<Rp;i++)
if(dis[i]<dis[k])
k=i;
printf("%.0lf\n",dis[k]);
}
return 0;
}模擬退火樣例2：
函式最值問題求解
題意：給出方程
分析：本題可以用經典的二分法求解，這種方法比較簡單，就不說了。主要來說模擬退火做法。之所以能夠二分，是因為F（x）的導數是單調的，所以只有一個極小值，所以可以直接二分。
程式碼1：
模擬退火大致步驟如下，類似樣例1的程式碼3,4：
1、隨機選取一個合適的控制條件T作為開始
2、隨機選取P個起始點，作為可行解
3、分別依據內能更新這P個可行解
4、減小控制條件T，直到終止條件
#include <iostream>
#include <string.h>
#include <stdlib.h>
#include <algorithm>
#include <stdio.h>
#include <time.h>
#include <math.h>
#define ITERS 10
#define T 100
#define eps 1e-8
#define delta 0.98
#define INF 1e99
using namespace std;
double x[ITERS];
int Judge(double x, double y)
{
if(fabs(x - y) < eps) return 0;
return x < y ? -1 : 1;
}
double Random()
{
double x = rand() * 1.0 / RAND_MAX;
if(rand() & 1) x *= -1;
return x;
}
double F(double x, double y)
{
return 6 * x * x * x * x * x * x * x + 8 * x * x * x * x * x * x + 7 * x * x * x + 5 * x * x - y * x;
}
void Init()
{
for(int i = 0; i < ITERS; i++)
x[i] = fabs(Random()) * 100;
}
double SA(double y)
{
double ans = INF;
double t = T;
while(t > eps)
{
for(int i = 0; i < ITERS; i++)
{
double tmp = F(x[i], y);
for(int j = 0; j < ITERS; j++)
{
double _x = x[i] + Random() * t;
if(Judge(_x, 0) >= 0 && Judge(_x, 100) <= 0)
{
double f = F(_x, y);
if(tmp > f)
x[i] = _x;
}
}
}
t *= delta;
}
for(int i = 0; i < ITERS; i++)
ans = min(ans, F(x[i], y));
return ans;
}
int main()
{
int t;
scanf("%d", &t);
while(t--)
{
double y;
scanf("%lf", &y);
srand(time(NULL));
Init();
printf("%.4lf\n", SA(y));
}
return 0;
}
#include<cstring>
#include<string>
#include<iostream>
#include<queue>
#include<cstdio>
#include<algorithm>
#include<map>
#include<cstdlib>
#include<cmath>
#include<vector>
//#pragma comment(linker, "/STACK:1024000000,1024000000");
using namespace std;
#define INF 0x3f3f3f3f
double y;
double getsum(double x)
{
return 6*pow(x,7)+8*pow(x,6)+7*pow(x,3)+5*x*x-y*x;
}
int main()
{
int T;
scanf("%d",&T);
while(T--)
{
scanf("%lf",&y);
double delte=0.98;
double t=100;
double x=0;
double ans=getsum(x);
while(t>1e-8)
{
int flag=1;
while(flag)
{
flag=0;
double temp=x+t;
if(temp>=0&&temp<=100&&getsum(temp)<ans&&fabs(ans-getsum(temp))>=1e-8)
{
x=temp;
ans=getsum(temp);
flag=1;
}
temp=x-t;
if(temp>=0&&temp<=100&&getsum(temp)<ans&&fabs(ans-getsum(temp))>=1e-8)
{
x=temp;
ans=getsum(temp);
flag=1;
}
}
t*=delte;
}
printf("%.4f\n",ans);
}
return 0;
}還有其他樣例，但是都沒搞懂怎麼寫模擬退火的過程，感覺沒有一個比較明確的方法，所以待明白之後再繼續深究，我只是寫了一個樣例1的一個程式，過了，我也不知為什麼可以過的，和樣例1程式碼差不多，只是
z.x = s.x + dx[i] * t;
z.y = s.y + dy[i] * t; z.x = s.x + dx[i] ;
z.y = s.y + dy[i] ;參考連結：（加*比較重要的連結或者沒看完的連結）