[hdu 5017 Ellipsoid] 模拟退火

[hdu 5017 Ellipsoid] 模拟退火

1. 题目链接

[hdu 5017 Ellipsoid]

2. 题意描述

给定一个三维空间的椭球面方程，求椭球面上的点到原点（0，0，0）的最小距离。

3. 解题思路

可以发现，椭球面上到原点的距离，具有一个极大值点和一个极小值点。
用模拟退火的方法可以近似搜索到全局最小。
这里因为只有一个极小值点，所以这里也不需要以一定概率接受比当前更差的解了。可以说是简化的模拟退火了。
一定程度上来说，模拟退火理论上可以近似替代二分、三分了。

4. 实现代码

#include <cmath>
#include <vector>
#include <string>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <algorithm>

using namespace std;

double a, b, c, d, e, f;
const int PSZ = 10;
const int inf = 0x3f3f3f3f;
const long long infl = 0x3f3f3f3f3f3f3f3f;
const double pi = acos(-1.0);
const double eps = 1e-8;
const int dir_sz = 8;

struct Point{
    double x, y, d;
} pt[PSZ];

double get_dis(double x, double y) {
    double ret = inf;
    double a1 = c;
    double b1 = d * y + e * x;
    double c1 = a * x * x + b * y * y + f * x * y - 1;
    double delta = b1 * b1 - 4 * a1 * c1;
    if(delta < 0) return ret;
    static vector<double> zs(2);
    zs[0] = (-b1 + sqrt(delta)) / (2 * a1);
    zs[1] = (-b1 - sqrt(delta)) / (2 * a1);
    for(double& z : zs) {
        double r = a * x * x + b * y * y + c * z * z + d * y * z + e * x * z + f * x * y;
        if(fabs(r - 1) > eps) continue;
        //printf("[%f %f %f %f]\n", x, y, z, r);
        ret = min(ret, sqrt(x * x + y * y + z * z));
    }
    return ret;
}

int main() {
#ifdef __LOCAL_WONZY__
    freopen("input.txt", "r", stdin);
#endif
    while(cin >> a >> b >> c >> d >> e >> f) {
        srand((unsigned int) time(NULL));
        pt[0].x = pt[1].x =  0;
        pt[0].y = sqrt(1 / b); pt[1].y = -sqrt(1 / b);
        pt[2].x = sqrt(1 / a); pt[3].x = -sqrt(1 / a);
        pt[2].y = pt[3].y = 0;
        double T = 1e3, r = 0.98, res = inf;
        int cnt = 0;
        for(int i = 0; i < 4; ++i) {
            res = min(res, pt[i].d = get_dis(pt[i].x, pt[i].y));
        }
        while(T > eps) {
            for(int i = 0; i < 4; ++i) {
                for(int j = 0; j < dir_sz; ++j) {
                    double x = pt[i].x + cos(2 * pi * j / dir_sz) * T;
                    double y = pt[i].y + sin(2 * pi * j / dir_sz) * T;
                    // ax^2 + by^2 + fxy = 1;
                    if(fabs(a * x * x + b * y * y + f * x * y - 1) < eps) continue;
                    double d = get_dis(x, y);
                    if(fabs(d - inf) < eps) continue;
                    double rand_db = (double) rand() / RAND_MAX;
                    double dE = d - pt[i].d;
                    res = min(res, d);
                    //只有一个极小值点，所以不需要一定概率接受比当前更差的解
                    //if(dE <= 0 || exp(-dE / T) > rand_db) {
                    if(dE <= 0) {
                        pt[i].x = x, pt[i].y = y, pt[i].d = d;
                    }
                }
            }
            T *= r;
            ++ cnt;
        }
        //cout << cnt << endl;
        cout << fixed << setprecision(7) << res << endl;
    }
    return 0;
}

5. 参考链接

1. GuoJiaSheng，模拟退火，https://www.cnblogs.com/GuoJiaSheng/p/4192301.html