算法设计课程项目报告

问题描述及要求

问题描述

Capacitated Facility Location Problem，限容量设施分配问题，是一类优化问题。目的是求解出代价最小的设施分配情况，包括设施的开关状态和服务顾客的情况。

题目要求

使用至少两种算法去求解问题。
使用71个测例进行测试，输出相应的结果。

问题分析

输入的数据包括每个设施的开设费用、每个设施的最大容量、每个顾客的需求，以及每个设施对应每个顾客的服务代价。以上三个是最核心的信息，我们需要求解出最优方案，使得在能满足所有顾客的需求下，并且在不超出设施的最大容量条件下，开设设施与服务顾客的总代价最小。

根据上面的分析，可以把该优化问题（$n$个设施，$m$个顾客）表达成数学形式：
$$
min \sum^n_{i=1}\sum^m_{j=1}c_{ij}y_{ij}+\sum^n_{i=1}f_ix_i \
s.t. \space\space \sum^n_{i=1}y_{ij} = 1 \space\text{for all } j=1,\dots, m \
\sum^m_{j=1}d_iy_{ij} \le u_ix_i \space \text{for all } i =1,\dots,m \
y_{ij} \ge 0 \space \text{for all } i = 1,\dots,m \
x_i \in \lbrace0,1\rbrace \space \text{for all }i = 1,\dots,n \
y_{ij} \in \lbrace0,1\rbrace \space \text{for all }j = 1,\dots,m
$$
其中的变量解释如下：

$x_i$表示第$i$个设施的开关状态，$x_i=0$表示设施关闭，$x_i=1$表示设施开启。
$f_i$表示第$i$个设施开启的费用。
$d_{j}$表示第$j$个顾客的需求。
$c_{ij}$表示第$i$个设施服务第$j$个顾客所需要的费用。
$y_{ij}$表示$i$个设施服务第$j$个顾客需求的比例（本次实验中，每位顾客只能由一个设施服务，因此$y_{ij}$要么为0，要么为1）。

分析：

我们极小化的部分由两部分组成，第一部分是设施服务顾客的花费，第二部分是设施开设的花费。
第一个约束条件的意思是，需要满足所有顾客的需求。
第二个约束条件的意思是，每个设施服务顾客的需求不能超出该设施的容量。

然后我们需要给出的输出有：

$x_i$：每个设施对应的开关状态。
每个客户选择的设施。

算法设计

方法1——贪心算法

思路：使用贪心算法求解非常简单，对于优化问题，贪心算法普遍适用，其核心思想是：为了达到最小值或最大值，我们在每一步选择的时候都是最贪新的。对于这个具体问题，我们的贪心策略是：对于每一个顾客，我们选取设施$i$，要求设施$i$服务这个顾客的代价是最小的。当然，这里的前提是设施$i$有足够容量服务这个顾客。于是每个顾客都按照这种贪心策略选择，最后就能得到一个最小代价了。我们可以根据顾客的选择来确定设施的开关，有顾客选择的设施就开放，没有顾客选择的设施就关闭。

贪心算法过程如下：

for all customers j
	if (presentCapacity[i] + demand[j] <= capacity[i]) {
		choose i that min(assignmentCost[i][j])
		update presentCapacity[i] += demand[j]
	}
	
for all facility i
	if (presentCapacity[i] != 0)
		openState[i] = 1

方法2——模拟退火法

思路：模拟退火法是元启发式搜索的其中一种方法。模拟退火法是克服爬山法缺点的有效方法，所谓退火是冶金专家为了达到某些晶体结构重复将金属加入或冷却的过程，该过程的控制参数为温度T。模拟退火法的基本思想是，在系统朝着能量减少的趋势这样一个变化过程中，偶尔允许系统跳到能量较高的状态，以避开局部极小点，最终稳定到全局最小点。

算法步骤如下：

随机挑选出一单元，并给他一个随机的位移$k$，求出系统因此而产生的能量变化$\Delta E_k$。
若$\Delta E_k \leq0$，该位移可采纳，而变化后的系统状态可作为下次变化的起点；若$\Delta E_k \gt0$，位移后的状态可采纳的概率为

$$
P_k = \frac{1}{(1+e^{-\Delta E_k/T})}
$$

式中T为温度，然后从(0,1)区间均匀分布的随机数中挑选一个数R。若$R \lt P_k$，则将变化后的状态作为下次的起点；否则，将变化前的状态作为下次的起点。

转第（1）步继续执行，直至达到平衡状态为止。

模拟退火法的关键在于如何产生新解，这里我采用了几种邻域搜索方式，目的是提高解的多样性。这里的解是设施的开启与关闭状态（虽然也可以使用设施与顾客的分配关系来作为解，但效果不佳）。方法有三个：

将某个位置的状态置反（若该设施开放，则改为关闭；反之亦然）
选出两个位置，倒转它们之间的状态。
选出两个位置，交换它们。

每次产生新解后，就要确定设施与顾客的分配关系，这个也是这次我做优化的地方。先说最开始的做法，就是根据生成的设施开放状态，为每个用户匹配一个设施（不超出容量限制），而这个匹配过程是随机的，即容量不超出限制的设施都可能会被选中。而我优化的地方是，在这个匹配过程加入了贪心算法的思想，按照概率使用两种匹配方式：
- 第一种是，匹配服务代价最小的设施给该顾客。
- 第二种是，随机匹配。
优化的思路是，如果只是使用随机匹配，那么每个顾客很可能挑选一些服务代价很大的设施，这是没有必要的。但是如果全部使用贪心的思想，也会出现问题，如果每个顾客只选自己代价最小的，那么可能导致一个问题：如果一个顾客C在设施A的代价很低，而在其他设施的代价很高，而其余顾客在各个设施的代价都差不多，那么顾客C有可能会因为设施A首先被占据而不得不选择代价极高的其他设施，但其实最合理的分配方式是顾客C选择设施A而其他顾客选择其他的设施。

需要注意的是因为通过邻域搜索的方式产生的解可能不能使得所有顾客都被服务，因此需要在匹配的时候做一个验证，如果无法匹配则重新生成新解。

代码实现

基本数据结构

#ifndef FLP_H
#define FLP_H

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

class FLP {
public:
  FLP();
  void input(); // 读入数据
  void initialize(); // 初始化变量
  void generateSolution(); // 邻域搜索产生新解
  int getTotalCost(); // 计算总代价
  bool update(); // 纯贪心匹配
  bool update1(); // 贪心+随机匹配
  int findBestFacility(int customerIndex); // 贪心策略
  int findRandomOne(int customerIndex); // 随机策略
  void greedy(); // 贪心算法
  void SA(); // 模拟退火法
  void create_new(); // 产生新解
private:
  // 输入内容
  int facilityNum; // 工厂数量
  int customerNum; // 顾客数量
  vector<int> capacity; // 工厂容量
  vector<int> openCost; // 工厂开设费用
  vector<double> demand; // 顾客需求
  vector<vector<double>> assignmentCost; // [i,j]第i个工厂对第j个顾客的费用

  vector<int> presentCapacity; // 当前时刻工厂已分配容量

  // 输出内容
  vector<int> openState; // 工厂开放状态
  vector<int> totalCost; // 工厂总支出
  vector<int> facilityChosen; // 顾客选择的工厂
    
  // 约束条件
  vector<bool> isSatisfied; // 顾客是否满足
};

#endif

贪心算法

对于每一个顾客，在设施列表中找到服务代价最小的工厂（使用一个presentCapacity来记录设施已经被使用的容量）。然后更新该设施的presentCapacity和代价，将匹配记录加入到facilityChose，设置isSatisfied为真。最后根据匹配的结果，判断设施的开放与关闭状态。

void FLP::greedy() {
  // 初始化
  for (int i = 0; i < facilityNum; i++) {
    openState[i] = 0;
  }

  for (int j = 0; j < customerNum; j++) {
    int choose;
    int minCost = INT_MAX;
    for (int i = 0; i < facilityNum; i++) {
      // 容量足够
      if (presentCapacity[i] + demand[j] <= capacity[i]) {
        // 每个客户找出最小的花费的设施
        if (assignmentCost[i][j] < minCost) {
          minCost = assignmentCost[i][j];
          choose = i;
        }
      }
    }
    presentCapacity[choose] += demand[j];
    //cout << j << " choose " << choose << endl;
    totalCost[choose] += assignmentCost[choose][j];
    isSatisfied[j] = true;
    facilityChosen[j] = choose;
  }

  for (int i = 0; i < facilityNum; i++) {
    //assert(presentCapacity[i]);
    if (presentCapacity[i] != 0) {
      openState[i] = 1;
      totalCost[i] += openCost[i];
    }
    else {
      openState[i] = 0;
    }
  }

  cout << getTotalCost() << endl;

  for (int i = 0; i < facilityNum; i++) {
    cout << openState[i] << " ";
  }
  cout << endl;
  for (int j = 0; j < customerNum; j++) {
    cout << facilityChosen[j] << " ";
  }
  cout << endl;

  // 检查是否所有顾客都满足
  for (int j = 0; j < customerNum; j++) {
    assert(isSatisfied[j]);
  }
}

模拟退火法

多种邻域搜索产生新解，使用多种方法目的是保证解的多样性。避免生成大量的无效解，在模拟退火的时候，解初始化的时候也是使用随机生成的。每次产生按照一定比例来产生，包括了反转、部分倒序、交换，这是我认为比较好的邻域搜索方式。

generateSolution()：

void FLP::generateSolution() {
  std::fill(presentCapacity.begin(), presentCapacity.end(), 0);
  std::fill(totalCost.begin(), totalCost.end(), 0);
  double r1 = ((double)rand()) / (RAND_MAX + 1.0);
  double r2 = ((double)rand()) / (RAND_MAX + 1.0);
  int pos1 = (int)(facilityNum * r1); // 第一个交叉点位置
  int pos2 = (int)(facilityNum * r2); // 第二个交叉点位置
  
  double ratio = rand() % 100 / (double)101;
  if (ratio <= 0.33) {
    if (openState[pos1] == 1) {
      openState[pos1] = 0;
    }
    else {
      openState[pos1] = 1;
    }
  }
  else if (ratio <= 0.66) {
    if (pos1 > pos2) {
      int t = pos1;
      pos1 = pos2;
      pos2 = t;
    }
    vector<int> v = openState;
    int j = pos2;
    for (int i = pos1; i <= pos2; i++) {
      openState[i] = v[j];
      j--;
    }
  }
  else {
    int x = openState[pos1];
    openState[pos2] = openState[pos1];
    openState[pos1] = x; 
  }
  
  /*
  cout << "New Solution: " << endl;
  for (int i = 0; i < facilityNum; i++) {
    cout << openState[i] << " ";
  }
  cout << endl;
  */
}

匹配的过程结合了随机与贪心，当该解无法匹配出结果时，返回false表示该解不合法，然后重新生成。

贪心策略：

int FLP::findBestFacility(int customerIndex) {
  int temp = INT_MAX;
  int index = -1;
  for (int i = 0; i < facilityNum; i++) {
    // 工厂开放且容量不超出限制
    if (openState[i] == 1 && presentCapacity[i] + demand[customerIndex] <= capacity[i]) {
      if (assignmentCost[i][customerIndex] < temp) {
        temp = assignmentCost[i][customerIndex];
        index = i;
      }
    }
  }
  return index;
}

随机策略：

int FLP::findRandomOne(int customerIndex) {
  int temp = INT_MAX;
  vector<int> indices;
  for (int i = 0; i < facilityNum; i++) {
    // 工厂开放且容量不超出限制
    if (openState[i] == 1 && presentCapacity[i] + demand[customerIndex] <= capacity[i]) {
      indices.emplace_back(i);
    }
  }
  if (indices.size() == 0) {
    return -1;
  }
  double r1 = ((double)rand()) / (RAND_MAX + 1.0);
  int pos = (int)(indices.size() * r1);
  return indices[pos];
}

匹配：

bool FLP::update1() {
  // 根据开关状态更新工厂开设费用
  for (int i = 0; i < facilityNum; i++) {
    if (openState[i] == 1) {
      totalCost[i] = openCost[i];
    }
  }
  // 给每个用户分配工厂（贪心）
  for (int j = 0; j < customerNum; j++) {
    int index;
    int r = rand() % 100 / (double)101;
    if (r <= 0.4) {
      index = findBestFacility(j);
    }
    else {
      index = findRandomOne(j);
    }
    //cout << "index: " << index;
    if (index == -1) {
      return false;
    }
    //cout << "i: " << index << " choose j: " << j << endl;
    facilityChosen[j] = index;
    presentCapacity[index] += demand[j];
    isSatisfied[j] = true;
    // 加上顾客费用
    totalCost[index] += assignmentCost[index][j];
  }
  return true;
}

模拟退火的过程可以按照上面给出的伪代码直接编写，有固定的模式，我们需要额外使用几个变量来记录整个过程中比较好的解。

void FLP::SA() {
  double T;
  T = T0;
  int count = 0;
  double f1, f2, df;
  double r;
  int best = INT_MAX;
  vector<int> bestOpenState;
  vector<int> bestFacilityChosen;

  // 初始化工厂开放状态
  for (int i = 0; i < facilityNum; i++) {
    double ratio = rand() % 100 / (double)101;
    if (ratio < 0.5) {
      openState[i] = 0;
    }
    else {
      openState[i] = 1;
    }
  }
  vector<int> origin = openState;
  update();

  f1 = getTotalCost();
  while (T > T_end) {
    for (int i = 0; i < L; i++) {
      origin = openState;
      create_new();

      f2 = getTotalCost();
      df = f2 - f1;
      if (df >= 0) {
        r = ((double)rand()) / (RAND_MAX);
        // 保留原来的解
        if (exp(-df / T) <= r) {
          openState = origin;
        }
      }
    }
    T *= q;
    if (best > getTotalCost()) {
      best = getTotalCost();
      bestOpenState = openState;
      bestFacilityChosen = facilityChosen;
    }
    count++;
  }

  cout << best << endl;
  for (int i = 0; i < facilityNum; i++) {
    cout << bestOpenState[i] << " ";
  }
  cout << endl;
  for (int j = 0; j < customerNum; j++) {
    cout << bestFacilityChosen[j] << " ";
  }
  cout << endl;

  // 检查是否所有顾客都满足
  for (int j = 0; j < customerNum; j++) {
    assert(isSatisfied[j]);
  }
}

测试结果及分析

贪心算法

Result Table（每个测例对应的详细结果请看github，这里篇幅有限，不详细展示）

TestCase	Result	OpenState	Time(s)
1	9440	1 1 1 1 1 1 1 1 1 1	0
2	8126	1 1 1 1 1 1 1 1 1 1	0
3	10126	1 1 1 1 1 1 1 1 1 1	0
4	12126	1 1 1 1 1 1 1 1 1 1	0.001
5	9357	1 1 1 1 1 1 1 1 1 1	0
6	8061	1 1 1 1 1 1 1 1 1 1	0
7	10061	1 1 1 1 1 1 1 1 1 1	0
8	12061	1 1 1 1 1 1 1 1 1 1	0
9	9040	1 1 1 1 1 1 1 1 1 1	0
10	7726	1 1 1 1 1 1 1 1 1 1	0.001
11	9726	1 1 1 1 1 1 1 1 1 1	0
12	11726	1 1 1 1 1 1 1 1 1 1	0
13	12032	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
14	9180	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
15	13180	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
16	17180	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
17	12032	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
18	9180	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
19	13180	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
20	17180	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
21	12032	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
22	9180	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0.001
23	13190	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
24	17180	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
25	19191	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
26	16131	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
27	21531	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
28	26831	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
29	19305	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
30	16239	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0.001
31	21639	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
32	27039	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
33	19055	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
34	15989	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
35	21389	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
36	26789	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
37	19055	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
38	15989	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0
39	21389	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0.001
40	26789	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1	0.001
41	7226	1 1 1 1 1 1 1 1 1 1	0.001
42	9957	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0.001
43	12448	1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1	0
44	7585	1 1 1 1 1 1 1 1 1 1	0
45	9848	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
46	12639	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1	0.001
47	6634	1 1 1 1 1 1 1 1 1 1	0
48	9044	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
49	12420	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0	0
50	10062	1 1 1 1 1 1 1 1 1 1	0
51	11351	1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1	0
52	10364	1 1 1 1 1 1 1 1 1 1	0
53	12470	1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1	0.001
54	10351	1 1 1 1 1 1 1 1 1 1	0
55	11970	1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1	0.001
56	23882	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0.001
57	32882	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
58	53882	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
59	39121	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
60	23882	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
61	32882	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
62	53882	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
63	39121	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
64	23882	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
65	32882	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0.001
66	53882	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
67	39671	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
68	23882	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0.001
69	32882	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0.001
70	53882	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0
71	39121	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1	0

模拟退火法

Result Table（每个测例对应的详细结果请看github，这里篇幅有限，不详细展示）

TestCase	Cost	OpenState	Time(s)
1	9142	1 1 1 1 1 1 0 0 1 1	76.68
2	8104	1 1 1 1 1 1 1 0 1 1	91.207
3	9824	1 1 1 1 1 1 1 0 1 0	74.669
4	11261	1 1 1 1 1 0 1 0 1 0	74.648
5	9348	1 1 1 1 1 1 0 1 1 1	416.138
6	8061	1 1 1 1 1 1 1 1 1 1	227.434
7	9905	1 1 1 1 1 1 1 1 1 0	439.569
8	11705	1 1 1 1 1 1 1 1 1 0	308.728
9	8598	1 1 1 1 1 0 1 0 1 0	52.849
10	7617	1 1 1 1 1 0 1 1 1 0	77.064
11	9064	1 1 1 1 1 0 1 0 1 0	92.202
12	10301	1 1 1 0 1 0 1 0 1 0	99.958
13	8858	1 0 0 1 0 0 0 0 0 1 1 0 1 1 0 1 1 1 0 1	108.256
14	7323	0 1 0 1 0 0 1 0 0 0 1 0 0 1 1 1 1 1 0 1	76.881
15	9568	0 0 0 1 0 0 1 1 0 0 1 0 0 0 1 0 1 1 1 1	74.479
16	11181	0 0 0 1 0 0 1 1 0 0 0 0 0 0 1 0 1 1 1 1	76.161
17	8564	1 0 0 1 0 0 1 0 0 1 1 0 0 1 1 1 0 1 0 1	72.535
18	7210	0 0 0 1 0 0 1 0 0 1 1 0 0 1 1 1 1 1 0 1	71040
19	9237	0 1 0 0 0 0 0 0 0 1 1 0 0 1 1 1 1 1 0 1	68.601
20	11122	0 0 1 1 0 0 1 0 0 0 1 0 1 0 0 1 1 1 0 1	73.531
21	8400	1 1 0 1 0 0 0 0 0 1 1 0 0 1 0 1 0 1 0 0	71.645
22	7288	1 1 0 0 0 0 1 0 0 1 1 0 0 1 1 1 0 1 0 1	74.509
23	9029	0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 1 1 1 0 1	71.480
24	10388	1 1 0 0 0 0 1 0 0 0 1 0 0 1 0 1 0 1 0 0	61578
25	13740	1 0 1 0 0 1 0 0 0 0 1 1 0 0 1 0 1 1 1 0 1 1 1 0 1 0 0 0 0 0	265.393
26	12098	1 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1	270.044
27	15451	1 0 1 0 0 1 0 0 0 0 1 1 0 0 1 0 1 1 0 0 1 0 1 1 1 1 0 0 0 1	284.05
28	17008	1 0 1 0 0 1 0 1 0 1 1 1 0 0 1 0 0 0 0 1 1 0 0 0 1 0 0 0 0 0	244.474
29	15533	1 0 1 1 0 1 0 1 1 0 0 1 0 1 1 0 0 1 0 0 1 1 1 0 1 0 0 1 0 1	261.252
30	14211	1 1 1 1 0 1 0 1 1 0 0 1 0 1 1 1 1 1 0 0 1 0 1 0 1 1 1 1 0 1	260.377
31	17724	1 0 1 0 0 1 0 0 0 1 0 1 1 1 1 1 1 1 0 0 1 0 0 1 1 0 0 1 1 1	250.266
32	21229	1 1 1 0 0 1 0 1 1 1 1 1 0 1 1 0 0 1 1 0 1 0 1 0 1 0 0 1 0 1	253。642
33	13829	1 0 1 0 0 1 1 0 1 1 0 1 0 0 0 0 0 1 0 0 1 0 1 0 1 0 0 0 0 0	237.683
34	12672	1 1 1 1 1 1 1 0 1 0 0 1 0 1 1 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0	268.002
35	14662	1 0 0 1 1 1 1 0 1 0 0 1 0 0 1 1 0 1 0 0 1 0 0 0 1 0 0 0 1 0	264.448
36	17572	1 0 1 1 0 1 0 0 0 0 1 1 0 1 1 0 0 1 0 0 1 0 0 0 1 1 0 0 0 1	259.474
37	13253	1 1 1 0 0 1 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 1 0 0 0 0	250.279
38	10880	1 0 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0	234.847
39	14562	0 0 1 0 0 1 0 1 0 0 0 0 0 0 1 0 0 1 1 1 1 0 0 0 1 0 1 0 0 0	234.299
40	15364	1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 0 0 1	233.660
41	6810	1 0 1 1 0 1 1 0 0 0	120.979
42	5766	0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 1	88.349
43	5748	0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 1	102.054
44	7585	1 1 1 1 1 1 1 1 1 1	99.524
45	6453	1 0 0 1 1 1 0 0 1 1 0 0 0 1 0 1 0 0 0 0	107.339
46	6563	0 0 0 1 0 0 1 0 1 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 1 1	113.745
47	6634	1 1 1 1 1 1 1 1 1 1	139.258
48	6054	1 0 0 0 0 1 0 1 1 0 1 0 0 1 1 0 1 0 1 0	93.707
49	5776	0 0 1 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0	105.128
50	9454	1 0 1 1 0 1 1 1 0 1	154.679
51	8315	0 1 0 0 0 1 1 0 0 1 1 0 1 1 0 1 1 0 0 1	125.583
52	9588	0 1 0 1 1 1 1 1 1 1	118.560
53	9176	0 1 0 1 0 1 1 0 0 1 0 0 1 0 0 1 0 0 1 0	123.052
54	9702	1 0 1 1 1 1 1 0 1 1	218.124
55	8387	0 0 1 0 0 1 1 1 0 1 0 0 0 1 1 0 0 0 1 0	132.905
56	22364	1 0 1 1 0 1 1 1 1 1 1 1 1 0 0 1 1 0 0 1 0 1 1 1 1 0 0 1 1 1	565.715
57	27908	1 1 1 0 1 0 1 1 1 1 1 0 1 0 0 1 0 0 1 1 1 0 0 1 1 1 0 1 1 0	592.974
58	40216	0 0 0 1 0 1 0 0 0 1 1 1 1 0 1 1 0 1 1 1 0 1 0 0 1 1 0 1 1 1	717.276
59	30156	0 0 1 0 0 0 1 0 1 0 1 1 0 1 1 1 1 1 1 1 1 0 0 1 1 1 0 0 1 1	643.627
60	21513	1 0 1 1 0 0 0 1 1 0 1 0 1 0 1 0 0 1 0 1 0 1 1 1 1 1 1 1 1 1	358.082
61	26200	1 0 0 0 1 1 0 1 0 0 0 0 1 1 0 0 1 0 0 1 1 1 1 0 0 0 0 1 1 0	361.7
62	34908	0 0 0 1 1 1 1 1 1 0 1 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 0	342.133
63	27719	0 0 1 0 1 0 0 0 1 0 0 1 0 1 0 1 1 0 1 1 1 0 0 1 1 1 0 0 1 0	347.366
64	21422	1 0 0 1 1 1 0 1 1 1 1 1 1 0 1 1 0 0 1 1 0 1 1 0 1 1 0 1 1 0	338.509
65	26174	1 0 0 1 1 0 0 1 1 0 0 1 1 0 1 1 0 0 0 1 0 1 1 1 1 0 0 0 1 0	336.559
66	33251	0 0 0 1 1 1 0 1 1 0 0 1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 1 0	340.147
67	27503	1 1 0 0 1 0 1 0 0 0 0 0 0 1 1 0 1 1 1 0 0 0 0 1 1 1 0 0 1 0	342.237
68	21747	1 1 0 0 1 1 1 1 1 0 1 0 1 1 1 0 1 0 1 1 0 1 1 0 1 0 0 1 1 0	342.120
69	26129	1 0 1 0 1 0 0 0 1 0 1 0 1 1 1 0 1 0 1 1 0 0 0 0 0 1 0 1 1 0	338.716
70	35195	1 1 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 1 1 1 0 1 0 0 0 0 0 1 0 1	338.060
71	28630	1 1 1 0 1 0 1 0 1 0 1 0 1 0 0 1 0 1 1 1 0 0 0 1 1 1 0 1 1 0	337.380

分析

直接对比两种算法，从结果来说，模拟退火法的结果更接近最优解，而且提高很大。这也说明对于这种动态规划问题，贪心算法几乎不可能求得最优解。而从时间上看，模拟退火法的代价极高，如果需要提高精度，所需要的时间可能会更长，而贪心算法则基本上是瞬间完成的，因为它不涉及到迭代。

总结

第一次遇到设施选址问题（Facility Location Problem），要求使用两种解法，于是想到了这个学期讲了很多的贪心算法，因为求解最优问题使用贪心算法一般能得到一个不错的结果（不是最优的！）。然后第二种解法选择用模拟退火法。模拟退火法的形式很简单，难点在于如何产生新解，其中有一个很重要的原则就是保证解的多样性，所以要使用多种邻域搜索策略。为了优化算法，我加入了一些贪心策略，得出来的结果也比一般的模拟退火法要好。

当然，两种做法的结果与优化算法计算出来的最优值还有一定差距，但是我认为模拟退火法计算出来的结果也是比较满意的。通过这次实验，我学会了如何使用贪心算法和模拟退火法，也能够根据具体问题进行分析，对算法做出适当的优化。当然，为了高效测试，还学会了使用批处理，的确提高了编程的效率。

整个项目的代码以及测试资料都在Github中，请访问：https://github.com/leungyukshing/Facility-Location-Problem

参考资料

优化算法求解结果参考：http://www.di.unipi.it/optimize/Data/MEX.html?tdsourcetag=s_pcqq_aiomsg
数值优化——Facility Location Problems：https://scipbook.readthedocs.io/en/latest/flp.html
Facility Location Problem解释：https://en.wikipedia.org/wiki/Facility_location_problem