codeforces 2C Commentator problem (2) -- 模拟退火

前面介绍了普通的计算过程，感觉还是比较繁琐和考验细致。在网上查了一下，看见不少人说用模拟退火来做，于是学习了一下，做了一个模拟退火的版本。

关键点在于能量函数的选取，状态产生函数的选取，以及参数的选取。

• 首先，能量函数的选取一定要仔细保证正确（开始犯了这个错误，如果发现程序老是会收敛到固定但是不正确的范围，可以考虑检查一下能量函数）。
• 其次，对于状态产生函数，候选解的产生应当带有一点随机性，并尽可能地遍历解空间。带有随机性的好处是增加跳出局部解的概率，否则容易陷入按规则走下去的局部解，同时，增加跳出局部解的概率了也更可能遍历解空间。而遍历解空间是并非随意地产生候选解，该题中候选解可以说是无穷多的，随意地产生候选解将无法收敛，在该题中可以逐步地减少步长以达到收敛。
• 最后，参数地选取最明显的是对于接受概率。需要让能量和温度取得一个适当的平衡，否则，接受概率始终过大和始终过小等于确定地接受或拒绝。当然，同时还有降温速度和准确度地平衡等等。

当温度降到预定值后，可以又重复模拟过程，因为后来步长减小了可能还是难以跳入另一个局部解（如果两个局部解距离很远），增大求得全局最优解的概率。该题我在采用这个重复步骤前很容易陷入不够好的解（对应两个解的情况）。

其实，能用确定算法的时候还是用确定的算法较好，随机算法调试比较麻烦，建议把过程结果输入到一个文件，便于调参数的时候看到过程，查看接受概率步长等等是否合适。

#include<cstdio>
#include<iostream>
#include<cmath>
#include<iomanip>
#include<cstdlib>
#include<ctime>

using namespace std;

struct Point{
  double a;
  double b;
  int r;
};

Point pt[3];
const int xdir[4] = {0,1,-1,0};
const int ydir[4] = {1,0,0,-1};
bool ability = false;
int base = 10000;

double calE(Point p){
  double l[3];
  double minl;
  double e = 0.0;
  int i;
  const int co = 1;

  for(i=0;i<3;i++){
    l[i] = pow(p.a - pt[i].a, 2) + pow(p.b - pt[i].b, 2);
    l[i] = l[i]/pow(pt[i].r, 2);
  }
  minl = fmin(l[0],l[1]);
  minl = fmin(minl,l[2]);

  for(i=0;i<3;i++){
    e += (l[i]-minl)*co;
  }
  if(e<0.001){
    e += minl;
    base = minl;
    ability = true;
  }else{
    e += base;
  }

  return e;
}

double move(Point *nextP, const Point p, double step){
  int j,dir;
  double e,nextE;
  Point stepP;

  dir  = rand()%4;
  step *= ((float) rand()) / (float) RAND_MAX;
  stepP.a = p.a + xdir[dir]*step;
  stepP.b = p.b + ydir[dir]*step;

  e = calE(stepP);
  *nextP = stepP;
  nextE = e;

  return nextE;
}

int main(){
  int i,j,k;
  double t=100,alpha=0.998,step;
  Point p,stepP;
  Point bestp;
  double bestE,e,oe;
  bool saFlag;
  int heatTime = 20;

#if defined(DEBUG)||defined(TEST)
  freopen("data","r",stdin);
  freopen("result","w",stdout);
#endif

  cout.setf(ios::fixed);
  srand(time(NULL));

  for(i=0;i<3;i++){
    cin>>pt[i].a>>pt[i].b>>pt[i].r;
  }
  bestp.a = p.a = ( pt[0].a + pt[1].a + pt[2].a )/3;
  bestp.b = p.b = ( pt[0].b + pt[1].b + pt[2].b )/3;
  bestE = calE(p);

#ifdef DEBUG
  cout<<"after initial"<<endl;
#endif

  while(heatTime>0){
    t = 100;
    while(t>0.000001){
      step = t*50;
      saFlag = true;
      oe = calE(p);

      e = move(&stepP, p, step);
#ifdef DEBUG
      cout<<"-------------------------"<<endl;
      cout<<e<<"\t"<<oe<<"\t"<<bestE<<endl;
      cout<< setprecision(5)<<stepP.a<<"\t"<<stepP.b<<"\t"<<p.a<<"\t"<<p.b<<endl;
      cout<< setprecision(6)<<"t:\t"<<t<<endl;
#endif
      if(e<bestE){
        bestp = p  = stepP;
        bestE = e;
        saFlag = false;
      }else{
        if(e<oe){
          p  = stepP;
          saFlag = false;
        }else{
          if( exp((oe-e)/t) > ((float) rand()) / (float) RAND_MAX ){//oe<=e
            oe = e;
            p = stepP;
          }
        }
      }

      if(saFlag){
        t *= alpha;
      }
    }
    heatTime --;
  }

  if(ability){
    cout<< setprecision(5) <<bestp.a<<" "<<bestp.b<<endl;
  }

  return 0;
}