从两道题来讨论一下模拟退火算法

今天上网课听了“连橘猫都能学会的计算几何”，感觉自己连橘猫都不如了…

今天唯一听懂的新知识应该就是这模拟退火算法了。

模拟退火算法是模拟钢铁才淬炼的时候的退火的过程温度变化的算法。

钢铁在退火的时候，其中某一点的温度是在不断变化的，也就是反复横跳的。模拟退火算法模拟了这一过程，在模拟精度达到一定的时候，可以实现得到全局最优解。

主要从两道题来讲。

第一道是二维的最小覆盖圆，第二道是三维的最小覆盖球。

最小覆盖圆的问题有几种解法：

• 1、随机增量法
• 2、爬山算法
• 3、模拟退火算法

随机增量法在对点进行随机化之后，时间复杂度可以优化到O(N)附近。

然而我用的是模拟退火来做的，就很玄学。

第一题是HDU – 3932

求最小覆盖圆。

由于这题的精度只要求到0.1，因此我用我笨拙的模拟退火算法过了这题。事实上这道题我这样子写模拟退火是有问题的，运行效率并不高。

#include<iostream>
#include<random>
#include<math.h>
#include<time.h>
#include<algorithm>

using namespace std;

typedef long long ll;
const double eps = 1e-8;
const int maxn = 1e3 + 10;

const int inf = 0x3f3f3f3f;
const double start_T = 20000;//退火起始温度

#define sqr(x) ((x) * (x))

bool dcmp(double a, double b)
{
    return fabs(a - b) < eps;
}

int X, Y;//地图的范围
int n;//其他点的总个数

default_random_engine engine(static_cast<unsigned int>(time(0)));



class point //点或向量的类
{
public:
    double x, y;
    void read()
    {
        cin >> x >> y;
    }
    void print()
    {
        cout << x << y;
    }

    point(double mx, double my)
    {
        x = mx;
        y = my;
    }

    point()
    {
        x = 0; y = 0;
    }

    point operator+(const point& p) const
    {
        return { x + p.x, y + p.y };
    }
    point operator-(const point& p) const
    {
        return { x - p.x, y - p.y };
    }
    point operator*(double p) const
    {
        return { x * p,y * p };
    }
    point operator/(double p)const
    {
        return { x / p,y / p };
    }
    double operator*(const point& p)
    {
        return x * p.x + y * p.y;
    }
    double len()
    {
        return sqrt(x * x + y * y);
    }
    double dis(const point& p)const
    {
        return ((*this) - p).len();
    }
    double angle()
    {
        return atan2(y, x);
    }



}pt[maxn];

double getMax(point p)
{
    double res = 0;
    for (int i = 0; i < n; i++)
    {
        res = max(res, pt[i].dis(p));
    }
    return res;
}

double SA()//模拟退火
{
    uniform_real_distribution<double> random_real_x(0, X);
    uniform_real_distribution<double> random_real_y(0, Y);
    uniform_real_distribution<double> random_real_tmp(0.0, 1.0);
    double T = start_T, rate = 0.94;
    point p;
    point to;
    p.x = random_real_x(engine);
    p.y = random_real_y(engine);

    double fx, fy, ans = getMax(p);

    int dx[] = { 1, 1, 1, -1, -1, -1, 0, 0 },
        dy[] = { 0, 1, -1, 0, 1, -1, -1, 1 };

    while (T > eps)
    {
        double tmp = inf;
        for (int times = 1; times <= 20; times++)
        {
            for (int i = 0; i < 8; i++)
            {
                //随机把点进行偏移，并且满足温度越低，偏移量越小
                fx = p.x + (dx[i] / start_T) * T * random_real_x(engine);
                fy = p.y + (dy[i] / start_T) * T * random_real_y(engine);

                if (fx < 0) fx = 0;
                if (fy < 0) fy = 0;
                if (fx > X) fx = X;
                if (fy > Y) fy = Y;

                double d = getMax(point(fx, fy));
                if (d < tmp)
                {
                    tmp = d;
                    to = point(fx, fy);
                }
            }
        }

        T *= rate;
        if (tmp < eps) continue;
        if (tmp < ans || random_real_tmp(engine) < exp(ans - tmp) / T * start_T)
        {
            ans = tmp;
            p = to;
        }
    }
    printf("(%.1f,%.1f).\n", p.x, p.y);
    return ans;
}


int main()
{

    while (cin >> X >> Y >> n)
    {
        for (int i = 0; i < n; i++)
        {
            //scanf_s("%lf%lf", pt[i].x, pt[i].y);
            cin >> pt[i].x >> pt[i].y;
        }
        printf("%.1f\n", SA());

    }

}

其实完全没必要像我这样去模拟退火，只需要随机偏移量然后直接移动就可以了。

但是，这样子的做法也不是最优化的。还有更优化的算法，完全不需要用到随机数。也就是每次退火过程都遍历已知点，寻找距离当前圆心最远的那个点，然后向那个方向移动一段距离，接着重复上述操作。经过多次上述操作，最终得到的点就是所求的最小覆盖圆的圆心。

上面这个方法是我做第二题的时候被卡时间了，才百度出来的。

第二题是三维的最小覆盖球的问题，题号是POJ – 2069

用上面代码的方法去做的话，就会卡时间TLE，毕竟我经过多次实验发现，其实这样子漫无目的的随机移动，虽然经过多次移动，还是能到达最终结果，但是这样的话精度很不稳定，而且耗时也很不稳定，从一百毫秒到一秒八都有。究其原因，是在生成随机数的时候很费时间而且做了很多的无谓的运算。因此就有了前面讲的定向移动的办法。这样子虽然不是完全模拟了退火，但是也满足了温度高的时候移动距离大，温度低移动距离小的条件。并且省去了生成随机数的过程，运算速度高了很多。

思路大概是这样子，然后就是代码实现了。

上代码之前不得不吐槽一下POJ，居然不支持C++11标准，害得我把代码改成了1998的标准，这POJ真的很不友好！！！

上代码！

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<math.h>
#include<time.h>
using namespace std;

typedef long long ll;
const double eps = 1e-10;
const int maxn = 40;
#define sqr(x) ((x) * (x))



bool dcmp(double a, double b)
{
    return fabs(a - b) < eps;
}

//三维点
struct point
{
    double x,y,z;
    void read()
    {
        scanf("%lf%lf%lf",&x,&y,&z);
    }
    

};

point pt[maxn];
int n;


double getMax(point p)
{
    double ans=0.0;
    for(int i=0;i<n;++i)
    {
        ans = max(ans,sqrt(sqr(p.x-pt[i].x)+sqr(p.y-pt[i].y)+sqr(p.z-pt[i].z)));
    }
    return ans;
}

double dis(point a,point b)
{
    return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y)+sqr(a.z-b.z));
}

const int inf = 0x3f3f3f3f;
const double start_T = 100.0;//设置起始温度
double rd()
{
    
    return 100*rand()/double(RAND_MAX);
}

//模拟退火
double SA()
{
    
    int dx[] = { 1, 1, 1, -1, -1, -1, 0, 0 },
        dy[] = { 0, 1, -1, 0, 1, -1, -1, 1 },
        dz[] = { -1,0, -1, 1, 0, 1,1,-1};
    
    double T = start_T;
    double rate = 0.98;//退火速度
    point p;
    p.x = 0; p.y = 0;p.z = 0;
    point to;

    double fx,fy,fz;
    double ans = inf;
    
    int s=0;
    while(T>eps)
    {
        double tmp = inf;
        for(int i=0;i<n;++i)//找出离当前点最远的那个点，以此来确定移动的方向。
        {
            if(dis(p,pt[s])<dis(p,pt[i])) s = i;
        }
        double mt = dis(p,pt[s]);
        ans = min(ans,mt);
        //向最远点移动以实现缩小r的目的
        p.x+=(pt[s].x-p.x)/mt*T;
        p.y+=(pt[s].y-p.y)/mt*T;
        p.z+=(pt[s].z-p.z)/mt*T;
        T*=rate;
        

    }

    return ans;

}


int main()
{
    
    while(true)
    {
        srand((int)time(0));
        //cout<<rd()<<endl;
        cin>>n;
        
        if(n==0) return 0;
        for(int i=0;i<n;i++)
        {
            pt[i].read();
        }
        

        printf("%.5lf\n",SA());
        //cout<<"耗时"<<end-start<<endl;
    }
}

欢迎关注我的公众号“灯珑”，让我们一起了解更多的事物~