模拟退火(SA)算法求解Max-Minsum Dispersion Problem（附代码及详细注释）

#include<iostream>
#include<cstdlib>
#include<cmath>
#include<string>
#include<ctime>
const int MAX = 0x7fffffff;

const int N = 1000;  //最大的范围
const int M = 500;  //要选择的集合大小
const int K = 100;  //两点间距离的最大值为K（距离默认为1-K）
const int max_count = 10; //当前温度的最大迭代次数
const double T0 = 50000.0; //初始温度
const double T_end = 1e-8; //退火结束温度
const double q = 0.98; //退火系数

int* elements; //共计N个点
int** distance;  //距离矩阵

clock_t start_total, end_total;  //计时器，整个程序
clock_t start_delta, end_delta;  //计时器，直接计算delta的步骤

struct Solution //解
{
 int* s0; //未被选中的数
 int* s1; //被选中的数
 int* delta; //到其他s1中的数的距离和
 int smallestDelta; //最大的delta，及目标函数值
 int center; //核心数
}iniSolution, bestSolution,solution1;

//分配存储空间
void init_solution(Solution* s)
{
 s->s0 = new int[M];
 s->s1 = new int[N - M];
 s->delta = new int[M];
 s->smallestDelta = 0;
}

//撤销iniSolution, bestSolution,solution1所占存储空间
void dispose(Solution* s)
{
 delete[](s->s0);s->s0 = NULL;
 delete[](s->s1);s->s1 = NULL;
 delete[](s->delta);s->delta = NULL;
}

//深拷贝solution类型
void copy_solution(Solution* ini, Solution* obj)
{
 for (int i = 0; i < M; i++)
  obj->s1[i] = ini->s1[i];
 for (int i = 0; i < M; i++)
  obj->s0[i] = ini->s0[i];
 for (int i = 0; i < M; i++)
  obj->delta[i] = ini->delta[i];
 obj->smallestDelta = ini->smallestDelta;
 obj->center = ini->center;
}

//计算所有delta的值
void calculate_delta(Solution* s)
{
 for (int i = 0; i < M; i++)
 {
  s->delta[i] = 0;
  for (int j = 0; j < M; j++)
   s->delta[i] += distance[s->s1[i]][s->s1[j]];
 }
}

// 在所有delta中找出smallest delta以及对应的中心数
void calculate_sum(Solution* s)
{
 s->smallestDelta = s->delta[0];
 s->center = s->s1[0];
 for (int i = 0; i < M; i++)
 {
  if (s->delta[i] < s->smallestDelta)
  {
   s->smallestDelta = s->delta[i];
   s->center = s->s1[i];
  }
 }
}

//更新的方法算出delta的值(将s0[v]与s1[u]交换）
void update_delta(Solution* s, int u, int v, int deltav)
{
 for (int i = 0; i < M; i++)
 {
  if (i == u)  //其自身delta的改变
   s->delta[u] = deltav - distance[s->s0[v]][s->s1[u]];
  else   //其他delta需将与s1[u]的距离转换为与s0[v]的距离
   s->delta[i] = s->delta[i] - distance[s->s1[u]][s->s1[i]] + distance[s->s0[v]][s->s1[i]];
 }
}

void init()
{
 //为距离矩阵随机生成1-100的距离
 for (int i = 0;i < N;i++)
  for (int j = 0;j < i;j++) //因为距离矩阵是对称的
   distance[i][j] = distance[j][i] = rand() % K + 1; //距离为1-K
 for (int i = 0;i < N;i++)
  distance[i][i] = 0;

 //随机生成初始解
 for (int i = 0; i < N; i++)
  elements[i] = i;
 //洗牌算法打乱
 for (int i = 0; i < N; i++)
 {
  int index = rand() % (N - i) + i;
  if (index != i)
  {
   int temp = elements[i];
   elements[i] = elements[index];
   elements[index] = temp;
  }
 }

 //初始化，分配数组空间
 init_solution(&iniSolution);
 //前M个为s1，后面为s0
 for (int i = 0;i < M;i++)
  iniSolution.s1[i] = elements[i];
 for (int i = M, j = 0;i < N;i++, j++)
  iniSolution.s0[j] = elements[i];

 //计算delta
 start_delta = clock();
 calculate_delta(&iniSolution);
 end_delta = clock();
 //计算smallest_delta
 calculate_sum(&iniSolution);
 //bestSolution拷贝iniSolution
 init_solution(&bestSolution);
 copy_solution(&iniSolution, &bestSolution);
 dispose(&iniSolution);
 //for (int i = 0;i < M;i++) std::cout << bestSolution.s1[i] << std::endl;
}

void SA_search() //模拟退火算法Simulated Annealing
{
 srand((unsigned)time(NULL)); //初始化随机数种子
 double T = T0; //初始温度
 int count_total = 0; //记录降温次数
 while (T > T_end) // 当温度低于结束温度时，退火结束
 {
  for (int count = 0;count <= max_count;count++) //count记录当前温度迭代次数
  {
   int deltav = 0; //计算deltav
   //产生新解solution1
   init_solution(&solution1);
   copy_solution(&bestSolution, &solution1);
   double r1 = ((double)rand()) / (RAND_MAX + 1.0);
   double r2 = ((double)rand()) / (RAND_MAX + 1.0);
   int v = (int)((N - M) * r1); //s0中交换点的位置
   int u = (int)(M * r2); //s1中交换点的位置
   for (int u = 0;u < M;u++) //对选中的数（s1）进行循环
    deltav += distance[bestSolution.s0[v]][bestSolution.s1[u]];
   update_delta(&solution1, u, v, deltav);

   int temp = solution1.s0[v];
   solution1.s0[v] = solution1.s1[u];
   solution1.s1[u] = temp;

   calculate_sum(&solution1); //计算smallest_delta
   double f1, f2, df;
   f1 = bestSolution.smallestDelta;
   f2 = solution1.smallestDelta;
   df = f2 - f1;
   double r = ((double)rand()) / (RAND_MAX); //0-1之间的随机数，用来决定是否接受新解
   if (df >= 0)
    copy_solution(&solution1, &bestSolution);
   else if (r < exp(df / T)) //若随机数小于p，接受新解
    copy_solution(&solution1, &bestSolution);
   dispose(&solution1);
   count++;
  }
  T *= q; //降温
  count_total++;
  std::cout << "第" << count_total << "次降温， 当前温度：" << T << "，当前最优解：" << bestSolution.smallestDelta << std::endl;
 }
}

void print_info()
{
 std::cout << "max-min sum answer:" << bestSolution.smallestDelta << std::endl;
 std::cout << "one delta run time:" << (double)(end_delta - start_delta) / CLOCKS_PER_SEC << std::endl;
 std::cout << "total run time:" << (double)(end_total - start_total) / CLOCKS_PER_SEC << std::endl;
 std::cout << "模拟退火算法，初始温度T0=" << T0 << "，降温系数q=" << q << ",每个温度迭代" << max_count << "次" << std::endl;
}

int main()
{
 //初始化数组，分配空间
 elements = new int[N];
 distance = new int* [N];
 for (int i = 0;i < N;i++)
  distance[i] = new int[N];

 init();
 start_total = clock();
 SA_search(); //模拟退火算法进行搜索
 end_total = clock();
 //打印结果
 print_info();

 dispose(&bestSolution);
 delete[]elements;
 for (int i = 0;i < N;i++)
  delete[]distance[i];
 delete[]distance;
 system("pause");
 return 0;
}