hdu 5017 Ellipsoid 模擬退火演算法 西安網路賽
阿新 • • 發佈：2019-02-05
For each test contains one line. Describes the minimal distance. Answer will be considered as correct if their absolute error is less than 10-5.
比賽的時候想到的是用拉格朗日乘數法求解附加條件下的多元函式的極值點。可惜列了四個方程組，這個是帶交叉項的，橢球經過了轉軸變換，三個軸都不平行於座標軸，簡直解不出來，太複雜了。
一.爬山演算法 ( Hill Climbing )
下面是兩種不同的方法：
比賽的時候想到的是用拉格朗日乘數法求解附加條件下的多元函式的極值點。可惜列了四個方程組，這個是帶交叉項的，橢球經過了轉軸變換，三個軸都不平行於座標軸，簡直解不出來，太複雜了。
一.爬山演算法 ( Hill Climbing )
介紹模擬退火前，先介紹爬山演算法。爬山演算法是一種簡單的貪心搜尋演算法，該演算法每次從當前解的臨近解空間中選擇一個最優解作為當前解，直到達到一個區域性最優解。
爬山演算法實現很簡單，其主要缺點是會陷入區域性最優解，而不一定能搜尋到全域性最優解。如圖1所示：假設C點為當前解，爬山演算法搜尋到A點這個區域性最優解就會停止搜尋，因為在A點無論向那個方向小幅度移動都不能得到更優的解。
圖1
二. 模擬退火(SA,Simulated Annealing)思想
爬山法是完完全全的貪心法，每次都鼠目寸光的選擇一個當前最優解，因此只能搜尋到區域性的最優值。模擬退火其實也是一種貪心演算法，但是它的搜尋過程引入了隨機因素。模擬退火演算法以一定的概率來接受一個比當前解要差的解，因此有可能會跳出這個區域性的最優解，達到全域性的最優解。以圖1為例，模擬退火演算法在搜尋到區域性最優解A後，會以一定的概率接受到E的移動。也許經過幾次這樣的不是區域性最優的移動後會到達D點，於是就跳出了局部最大值A。
模擬退火（Simulated Annealing，簡稱SA）是一種通用概率演算法 ，用來在一個大的搜尋空間內找尋命題的最優解。
“模擬退火”演算法是源於對熱力學中退火過程的模擬，在某一給定初溫下，通過緩慢下降溫度引數，使演算法能夠在多項式時間內給出一個近似最優解。退火與冶金學上的‘退火’相似，而與冶金學的淬火有很大區別，前者是溫度緩慢下降，後者是溫度迅速下降。
原理編輯
“模擬退火”的原理也和金屬退火的原理近似：我們將熱力學的理論套用到統計學上，將搜尋空間內每一點想像成空氣內的分子；分子的能量，就是它本身的動能；而搜尋空間內的每一點，也像空氣分子一樣帶有“能量”，以表示該點對命題的合適程度。演算法先以搜尋空間內一個任意點作起始：每一步先選擇一個“鄰居”，然後再計算從現有位置到達“鄰居”的方法
演算法對應動態演示圖： 模擬退火演算法新解的產生和接受可分為如下四個步驟： 第一步是由一個產生函式從當前解產生一個位於解空間的新解；為便於後續的計算和接受，減少演算法耗時，通常選擇由當前新解經過簡單地變換即可產生新解的方法，如對構成新解的全部或部分元素進行置換、互換等，注意到產生新解的變換方法決定了當前新解的鄰域結構，因而對冷卻進度表的選取有一定的影響。 第二步是計算與新解所對應的目標函式差。因為目標函式差僅由變換部分產生，所以目標函式差的計算最好按增量計算。事實表明，對大多數應用而言，這是計算目標函式差的最快方法。 第三步是判斷新解是否被接受,判斷的依據是一個接受準則，最常用的接受準則是Metropolis準則: 若Δt′<0則接受S′作為新的當前解S，否則以概率exp(-Δt′/T)接受S′作為新的當前解S。 第四步是當新解被確定接受時，用新解代替當前解，這隻需將當前解中對應於產生新解時的變換部分予以實現，同時修正目標函式值即可。此時，當前解實現了一次迭代。可在此基礎上開始下一輪試驗。而當新解被判定為捨棄時，則在原當前解的基礎上繼續下一輪試驗。 模擬退火演算法與初始值無關，演算法求得的解與初始解狀態S(是演算法迭代的起點)無關；模擬退火演算法具有漸近收斂性，已在理論上被證明是一種以概率l 收斂於全域性最優解的全域性優化演算法；模擬退火演算法具有並行性。下面是兩種不同的方法：
//78MS 296K
#include <cstdio>
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std;
double a, b, c, d, e, f;
const double INF = 1LL<<60;
const double EPS = 1e-8;
int dx[] = {-1, -1, -1, 0, 0, 1, 1, 1};
int dy[] = {-1, 0, 1, -1, 1, -1, 0, 1};
inline double dis(double x, double y, double z){
return sqrt(x*x + y*y + z*z);
}
inline double getz(double x, double y){
double ta = c, tb = d*y+e*x, tc = a*x*x+b*y*y+f*x*y-1;
double delta=tb*tb-4*ta*tc;
if(delta<0) return INF;
double z1=(-tb+sqrt(delta))/(2*ta);
double z2=(-tb-sqrt(delta))/(2*ta);
return z1*z1<z2*z2?z1:z2;
}
void solve(){
double x=0,y=0,z=getz(x,y),delta=0.8;
while(delta>EPS){
for(int i=0;i<8;i++){
double tx=x+dx[i]*delta,
ty=y+dy[i]*delta,
tz=getz(tx,ty);
if(tz>INF/10) continue;
if(dis(tx,ty,tz)-dis(x,y,z)<0)
x=tx,y=ty,z=tz;
}
delta *= 0.99;
}
printf("%.7f\n",dis(x,y,z));
}
int main(){
//srand((unsigned)time(NULL));
while(~scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&e,&f)){
solve();
}
return 0;
}
//78MS 296K
#include <cstdio>
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std;
double a, b, c, d, e, f;
const double INF = 1LL<<60;
const double EPS = 1e-8;
const int dx[] = {-1, -1, -1, 0, 0, 1, 1, 1};
const int dy[] = {-1, 0, 1, -1, 1, -1, 0, 1};
//int dx[] = {1, 0, -1, 0};
//int dy[] = {0, 1, 0, -1};
inline double dis(double x, double y, double z){
return sqrt(x*x + y*y + z*z);
}
inline double getz(double x, double y){
double ta = c, tb = d*y+e*x, tc = a*x*x+b*y*y+f*x*y-1;
double delta=tb*tb-4*ta*tc;
if(delta<0) return INF;
double z1=(-tb+sqrt(delta))/(2*ta);
double z2=(-tb-sqrt(delta))/(2*ta);
return z1*z1<z2*z2?z1:z2;
}
void solve(){
double x=0, y=0, z=getz(x,y), delta=0.8;
while(delta > EPS){
double nowd = dis(x, y, z);
double nxtx = x, nxty = y, nxtz = z;
for(int i=0;i<8;i++){
double tx=x+dx[i]*delta,
ty=y+dy[i]*delta,
tz=getz(tx,ty);
if(tz>INF/10) continue;
if(dis(tx, ty, tz) - nowd < 0){
nowd = dis(tx, ty, tz);
nxtx = tx, nxty = ty, nxtz = tz;
}
}
x = nxtx, y = nxty, z = nxtz;
delta *= 0.99;
}
printf("%.7f\n",dis(x,y,z));
}
int main(){
//srand((unsigned)time(NULL));
while(~scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&e,&f)){
solve();
}
return 0;
}