hdu 5017 Ellipsoid 模拟退火算法 西安网络赛

Given a 3-dimension ellipsoid(椭球面)

your task is to find the minimal distance between the original point (0,0,0) and points on the ellipsoid. The distance between two points (x₁,y₁,z₁) and (x₂,y₂,z₂) is defined as

 

Input

There are multiple test cases. Please process till EOF.

For each testcase, one line contains 6 real number a,b,c(0 < a,b,c,< 1),d,e,f(0 ≤ d,e,f < 1), as described above.It is guaranteed that the input data forms a ellipsoid. All numbers are fit in double.

 

Output

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 )

         介绍模拟退火前，先介绍爬山算法。爬山算法是一种简单的贪心搜索算法，该算法每次从当前解的临近解空间中选择一个最优解作为当前解，直到达到一个局部最优解。

         爬山算法实现很简单，其主要缺点是会陷入局部最优解，而不一定能搜索到全局最优解。如图1所示：假设C点为当前解，爬山算法搜索到A点这个局部最优解就会停止搜索，因为在A点无论向那个方向小幅度移动都不能得到更优的解。

图1

 

 

二. 模拟退火(SA,Simulated Annealing)思想

         爬山法是完完全全的贪心法，每次都鼠目寸光的选择一个当前最优解，因此只能搜索到局部的最优值。模拟退火其实也是一种贪心算法，但是它的搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。以图1为例，模拟退火算法在搜索到局部最优解A后，会以一定的概率接受到E的移动。也许经过几次这样的不是局部最优的移动后会到达D点，于是就跳出了局部最大值A。

模拟退火（Simulated Annealing，简称SA）是一种通用概率算法，用来在一个大的搜寻空间内找寻命题的最优解。

“模拟退火”算法是源于对热力学中退火过程的模拟，在某一给定初温下，通过缓慢下降温度参数，使算法能够在多项式时间内给出一个近似最优解。退火与冶金学上的‘退火’相似，而与冶金学的淬火有很大区别，前者是温度缓慢下降，后者是温度迅速下降。

原理编辑

“模拟退火”的原理也和金属退火的原理近似：我们将热力学的理论套用到统计学上，将搜寻空间内每一点想像成空气内的分子；分子的能量，就是它本身的动能；而搜寻空间内的每一点，也像空气分子一样带有“能量”，以表示该点对命题的合适程度。算法先以搜寻空间内一个任意点作起始：每一步先选择一个“邻居”，然后再计算从现有位置到达“邻居”的概率。

模拟退火算法可以分解为解空间、目标函数和初始解三部分。

模拟退火的基本思想:

(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)， 每个T值的迭代次数L

(2) 对k=1，……，L做第(3)至第6步：

(3) 产生新解S′

(4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数

(5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/(KT))接受S′作为新的当前解(k为波尔兹曼常数）.

(6) 如果满足终止条件则输出当前解作为最优解，结束程序。

终止条件通常取为连续若干个新解都没有被接受时终止算法。

(7) T逐渐减少，且T->0，然后转第2步。

方法

算法对应动态演示图：

模拟退火算法新解的产生和接受可分为如下四个步骤：

第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。

第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。

第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若Δt′<0则接受S′作为新的当前解S，否则以概率exp(-Δt′/T)接受S′作为新的当前解S。

第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。

模拟退火算法与初始值无关，算法求得的解与初始解状态S(是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。

下面是两种不同的方法：

//78MS296K#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;}

//78MS296K#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;}