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

Part 1 Max-Minsum Dispersion Problem

先来看一个小故事，转自(链接：http://blog.csdn.net/fudan_abc/article/details/2052642)，假如老板要你解决一个问题，你绞尽脑汁还是想不出来，叫天天不应，叫地地不灵，这时你走进老板办公室，可以采取3种策略：

(1)一副倒霉像，神情ws，可怜巴巴的说：老板，我没做出来，我想我是太蠢了。。。
boss：蠢材！滚!
(失败。。。)

(2)雄赳赳气昂昂跨进老板办公室，大吼一声：小样，你丫给我问题根本就无解，害我白想这么些天，我靠！
boss：我才靠，自己做不出来就说这个问题无解，要是人人都这样混，我这老板还当个屁啊，滚！
(做不出来还如此气愤，不仅失败，而且欠扁。。。)

(3)从容不迫的说：老板，我做不出来，但是，我敢肯定，那些大牛们也照样做不出来
boss：原来是这样，那也难为你了。

把决策问题可以按照难易程度分为几类：
(1)P问题：可以在多项式( polynomial )时间内解决的问题，称为P问题。
(2)NP问题：对于一类问题，我们可能没有一个已知的快速的方法得到问题的答案，但给定一个解，我们可以在P时间内检查他正确与否的决策问题，成为NP( Non-deterministic polynomial)问题。
(4)NP-hard问题：用一句话概括他们的特征就是“at least as hard as the hardest problems in NP Problem”。
【多项式级时间复杂度：O(1),O(log(n)),O(n^a)等。因为规模n出现在底数的位置。】

1.1 Max-Minsum Dispersion Problem

Max-Minsum DP就是一个典型的NP-Hard问题，属Equity-based dispersion problems，试图从较大的集合中选择一组元素时解决公平与效率的平衡，应用广泛然而计算难度较大。
简单来讲，就是要从一个集合中选择一个子集合，使得子集合中某个所选元素到其他所选元素之间距离的最小和最大化。

举个例子，假如说你一天要完成5个任务，现在有10项任务供你选择，煮饭30分钟，做菜30分钟，衣服机洗30分钟，做作业20分钟，烧开水10分钟，吃饭15分钟等，你是个聪明的无聊人，知道有些活儿可以一起干以节省时间，又想让这5项任务占用你更多的时间以打发无聊，这便是个简单的Max-Minsum DP。

再举个更贴近实际的，对于网页排名问题，即是标识网页等级或重要性。最早的搜索引擎采用的是分类目录的方法，即通过人工对网页进行分类并整理出高质量网站。随着网页数目的急剧增大，这种方法显然无法实现，Larry Page和Sergey Brin受学术界对学术论文重要性的评估方法（论文引用次数）的启发，提出了PageRank算法，如果一个网页被很多其它网页链接到，说明这个网页很重要，它的PageRank值也会相应较高，如果一个PageRank值很高的网页链接到另外某个网页，那么那个网页的PageRank值也会相应地提高。所以说找重要论文、相关论文就不免涉及到Max-Minsum DP。

更多的应用如下图：

1.2 Max-Minsum DP的数学描述

考虑一个含有n个元素的集合N={1,2,...,n}，每个元素包含着r个属性，我们可以将一个元素用向量表示。问题在于选择N的一个子集M，|M|为固定正整数m（m

这个距离有多种算法，如欧几里得距离，曼哈顿距离等。在这里我们使用最为常用的欧几里得距离

问题可以表达为：

Part 2 模拟退火算法(SA)再回顾

在之前的推文【算法进阶】用模拟退火(SA, Simulated Annealing)算法解决旅行商问题中，已经对模拟退火算法有了详细介绍并给出了伪代码及实例。在这里，我们简要复习，详细参见以上推文。

2.1 SA算法介绍

模拟退火算法的基础是metropolis算法。metropolis算法又称为metropolis抽样，其核心思想是：当能量增加的时候以一定概率接纳，而非一味拒绝。
所以，当Y(i+1)>Y(i)，则无条件接受；
   当Y(i+1)  以一定概率接受一个比当前解较差的解，从而在一定程度上避免陷入局部最优。
  然而应当如何计算这个概率呢？根据热力学的原理，在温度为T时，出现能量差为dE的降温的概率为P(dE)， 表示为：

Part 3 具体算法介绍

• 通过邻域动作产生新的集合M，产生新解及当前解，计算即smallestDelta。
• 若当前解的smallestDelta小于最优解的smallestDelta，则更新最优解为当前解，否则以模拟退火的那个概率接受当前解，然后降温。
• 重复之前步骤，直到满足退出条件。

现在拿一个小算例来操作一下：
现有一点集N={(0,1),(1,2),(3,4),(4,5),(6,6),(8,7)}，我们要从中选出m个点构成点集M，就取m=3吧，目标函数是我们挑选的这3个点中的某一点到其余2个点的距离之和的最小值，而问题在于找到使目标函数值最大的那3个点。

3.1 初始解生成

就小算例而言，我们就随机选3个点，你不妨可以掷骰子，我掷的是5，2，1，那我们就取(6,6), (1,2), (0,1)这三个点，不妨将这三个点重新标记为,
以为中心点，则;
以为中心点，则;
以为中心点，则;
不难看出，smallestDelta为Δ2，故初始解也是当前最优解即为M={(6,6),(1,2),(0,1)}，对应为。
  而就本问题而言，对于初始解，我们亦是随机产生，距离矩阵利用洗牌算法随机生成1-100的距离，随机选择m个元素构成s1，未被选中的即为s0，为了识别M和N\M，我们利用n维向量∆= (∆1 , ∆2 , . . . , ∆n )，其中，则，对应的集合M即为初始解，也作为最优解。

3.2 邻域动作

采用exchange算子：从被选择的元素的集合中随机选择元素u，即u∈M，从不被选择的元素的集合中随机选择元素v，即v∈N\M，交换u, v。拿上文小算例N={(0,1),(1,2),(3,4),(4,5),(6,6),(8,7)}举个例子，从中随机选择，即∈M，从中随机选择，即∈N\M，交换，此时得到新解，以三点分别为中心点，故不变得到
以为中心点，则;
以为中心点，则;
以为中心点，则;
不难看出，smallestDelta为Δ2，故最优解更新，变为M={(4,5),(1,2),(0,1)}，对应为。

3.3 去重优化

对于本问题，给定邻域解和对应向量(∆1 , ∆2 , . . . , ∆n )，目标值可以在O(M)时间内计算，此外，若是两个元素u∈M，v∈N\M交换，则向量∆= (∆1 , ∆2 , . . . , ∆n )可以在O(N)时间内快速更新，具体可表示为下图：

为了通俗易懂，接着拿上文小算例N={(0,1),(1,2),(3,4),(4,5),(6,6),(8,7)}举例，比较3.1及3.2计算Δ过程不难看出，对于未改变的点，即以为中心点、以为中心点时，对应的Δ计算过程只改变了一半，这部分就是我们可以优化的部分，因为另外一半我们就不用再重复计算了，

；
；
当数据越来越庞大之后，这部分优化带来的效益就会体现得更加明显，时间复杂度大幅减少。
而对于改变的点，
， 基本上就是重算。

Part 4 代码分享

算例为随机生成，具体实现如下：

#include
#include
#include
#include
#include
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;
}

结果如图：