一尘不染

优化的TSP算法

algorithm

我感兴趣的方法来改善或者想出了能够解决算法旅行商问题有关n = 100 to 200的城市。

我提供的Wikipedia链接列出了各种优化,但是这样做的水平很高,我不知道如何在代码中实际实现它们。

那里有工业强度求解器,例如Concorde,但是对于我想要的来说太复杂了,而泛滥对TSP的搜索的经典解决方案都提供了随机算法或经典的回溯或动态编程算法,这些算法只能用于20个城市。

那么,有谁知道如何实现一个至少在100个城市的合理时间(几秒钟)内可以运行的简单TSP求解器(简单地说,一个实现不需要100-200行代码)?我只对确切的解决方案感兴趣。

您可能会假设输入将是随机生成的,所以我不在乎专门针对破坏特定算法的输入。


阅读 357

收藏
2020-07-28

共1个答案

一尘不染

200行且无库是一个严格的约束。高级求解器使用分支定界和Held–Karp松弛约束,我不确定哪怕是最基本的版本也能容纳200条法线。不过,这是一个大纲。

举行卡普

将TSP编写为整数程序的一种方法如下(Dantzig,Fulkerson,Johnson)。对于所有边e,常数w
e表示边e的长度,如果边e在巡视线上,则变量x e为1,否则为0。对于所有顶点S的子集,∂(S)表示连接S中的顶点和非S中的顶点的边。

最小化总和e w e x e
服从
1.对所有顶点v,总和e在∂({v})中 x e = 2
2.对于所有非空适当顶点子集S,总和e在∂(S)中 X ë ≥2
3.所有边缘于E E,X ë在{0,1}

条件1确保边缘集是路线的集合。条件2确保只有一个。(否则,让S为其中一个巡回线访问的一组顶点。)通过进行此更改,可以使Held–Karp松弛。

3.对于所有的边缘E在E,X ë在{0,1}
3.在E,0≤X所有边缘ë ë ≤1

Held–Karp是一个线性程序,但是它的约束数量是指数的。解决该问题的一种方法是引入Lagrange乘法器,然后进行次梯度优化。归结为一个循环,该循环计算最小生成树,然后更新一些向量,但其中涉及到细节。除了“
Held–Karp”和“次梯度(下降)”之外,“ 1-tree”是另一个有用的搜索词。

(一种较慢的替代方法是编写一个LP解算器,并引入子行程约束,因为先前的最优方法违反了这些约束。这意味着编写LP解算器和最小切割程序,这也是更多代码,但可能会更好地扩展到更特殊的TSP约束。)

分支定界

所谓“部分解决方案”,是指变量的部分分配为0或1,其中分配为1的边肯定在游览中,分配为0的边肯定是出去。利用这些附带条件评估Held–Karp可以使最佳行程的下界符合已做出的决定(扩展)。

分支定界法维护一组局部解决方案,其中至少一个扩展到最佳解决方案。具有最佳优先回溯的一种变体,深度优先搜索的伪代码如下。

let h be an empty minheap of partial solutions, ordered by Held–Karp value
let bestsolsofar = null
let cursol be the partial solution with no variables assigned
loop
    while cursol is not a complete solution and cursol's H–K value is at least as good as the value of bestsolsofar
        choose a branching variable v
        let sol0 be cursol union {v -> 0}
        let sol1 be cursol union {v -> 1}
        evaluate sol0 and sol1
        let cursol be the better of the two; put the other in h
    end while
    if cursol is better than bestsolsofar then
        let bestsolsofar = cursol
        delete all heap nodes worse than cursol
    end if
    if h is empty then stop; we've found the optimal solution
    pop the minimum element of h and store it in cursol
end loop

分支定界的想法是,存在一个部分解决方案的搜索树。解决Held–Karp的要点是LP的值最多为最佳行程的OPT长度,但推测至少为3/4
OPT(实际上,通常更接近OPT)。

我遗漏的伪代码中的一个细节是如何选择分支变量。通常,目标是首先做出“艰难”的决定,因此,修复其值已经接近0或1的变量可能是不明智的。一种选择是选择最接近0.5的值,但还有很多其他的选择。

编辑

Java实现。198条非空白,非注释行。我忘记了1树不能将变量分配给1,所以我通过查找1树的度大于2的顶点进行分支并依次删除每个边。该程序接受TSPLIB实例的EUC_2D形式,例如,eil51.tspeil76.tspeil101.tsplin105.tsphttp://www2.iwr.uni-
heidelberg.de/groups/comopt/software/TSPLIB95/tsp/。

// simple exact TSP solver based on branch-and-bound/Held--Karp
import java.io.*;
import java.util.*;
import java.util.regex.*;

public class TSP {
  // number of cities
  private int n;
  // city locations
  private double[] x;
  private double[] y;
  // cost matrix
  private double[][] cost;
  // matrix of adjusted costs
  private double[][] costWithPi;
  Node bestNode = new Node();

  public static void main(String[] args) throws IOException {
    // read the input in TSPLIB format
    // assume TYPE: TSP, EDGE_WEIGHT_TYPE: EUC_2D
    // no error checking
    TSP tsp = new TSP();
    tsp.readInput(new InputStreamReader(System.in));
    tsp.solve();
  }

  public void readInput(Reader r) throws IOException {
    BufferedReader in = new BufferedReader(r);
    Pattern specification = Pattern.compile("\\s*([A-Z_]+)\\s*(:\\s*([0-9]+))?\\s*");
    Pattern data = Pattern.compile("\\s*([0-9]+)\\s+([-+.0-9Ee]+)\\s+([-+.0-9Ee]+)\\s*");
    String line;
    while ((line = in.readLine()) != null) {
      Matcher m = specification.matcher(line);
      if (!m.matches()) continue;
      String keyword = m.group(1);
      if (keyword.equals("DIMENSION")) {
        n = Integer.parseInt(m.group(3));
        cost = new double[n][n];
      } else if (keyword.equals("NODE_COORD_SECTION")) {
        x = new double[n];
        y = new double[n];
        for (int k = 0; k < n; k++) {
          line = in.readLine();
          m = data.matcher(line);
          m.matches();
          int i = Integer.parseInt(m.group(1)) - 1;
          x[i] = Double.parseDouble(m.group(2));
          y[i] = Double.parseDouble(m.group(3));
        }
        // TSPLIB distances are rounded to the nearest integer to avoid the sum of square roots problem
        for (int i = 0; i < n; i++) {
          for (int j = 0; j < n; j++) {
            double dx = x[i] - x[j];
            double dy = y[i] - y[j];
            cost[i][j] = Math.rint(Math.sqrt(dx * dx + dy * dy));
          }
        }
      }
    }
  }

  public void solve() {
    bestNode.lowerBound = Double.MAX_VALUE;
    Node currentNode = new Node();
    currentNode.excluded = new boolean[n][n];
    costWithPi = new double[n][n];
    computeHeldKarp(currentNode);
    PriorityQueue<Node> pq = new PriorityQueue<Node>(11, new NodeComparator());
    do {
      do {
        boolean isTour = true;
        int i = -1;
        for (int j = 0; j < n; j++) {
          if (currentNode.degree[j] > 2 && (i < 0 || currentNode.degree[j] < currentNode.degree[i])) i = j;
        }
        if (i < 0) {
          if (currentNode.lowerBound < bestNode.lowerBound) {
            bestNode = currentNode;
            System.err.printf("%.0f", bestNode.lowerBound);
          }
          break;
        }
        System.err.printf(".");
        PriorityQueue<Node> children = new PriorityQueue<Node>(11, new NodeComparator());
        children.add(exclude(currentNode, i, currentNode.parent[i]));
        for (int j = 0; j < n; j++) {
          if (currentNode.parent[j] == i) children.add(exclude(currentNode, i, j));
        }
        currentNode = children.poll();
        pq.addAll(children);
      } while (currentNode.lowerBound < bestNode.lowerBound);
      System.err.printf("%n");
      currentNode = pq.poll();
    } while (currentNode != null && currentNode.lowerBound < bestNode.lowerBound);
    // output suitable for gnuplot
    // set style data vector
    System.out.printf("# %.0f%n", bestNode.lowerBound);
    int j = 0;
    do {
      int i = bestNode.parent[j];
      System.out.printf("%f\t%f\t%f\t%f%n", x[j], y[j], x[i] - x[j], y[i] - y[j]);
      j = i;
    } while (j != 0);
  }

  private Node exclude(Node node, int i, int j) {
    Node child = new Node();
    child.excluded = node.excluded.clone();
    child.excluded[i] = node.excluded[i].clone();
    child.excluded[j] = node.excluded[j].clone();
    child.excluded[i][j] = true;
    child.excluded[j][i] = true;
    computeHeldKarp(child);
    return child;
  }

  private void computeHeldKarp(Node node) {
    node.pi = new double[n];
    node.lowerBound = Double.MIN_VALUE;
    node.degree = new int[n];
    node.parent = new int[n];
    double lambda = 0.1;
    while (lambda > 1e-06) {
      double previousLowerBound = node.lowerBound;
      computeOneTree(node);
      if (!(node.lowerBound < bestNode.lowerBound)) return;
      if (!(node.lowerBound < previousLowerBound)) lambda *= 0.9;
      int denom = 0;
      for (int i = 1; i < n; i++) {
        int d = node.degree[i] - 2;
        denom += d * d;
      }
      if (denom == 0) return;
      double t = lambda * node.lowerBound / denom;
      for (int i = 1; i < n; i++) node.pi[i] += t * (node.degree[i] - 2);
    }
  }

  private void computeOneTree(Node node) {
    // compute adjusted costs
    node.lowerBound = 0.0;
    Arrays.fill(node.degree, 0);
    for (int i = 0; i < n; i++) {
      for (int j = 0; j < n; j++) costWithPi[i][j] = node.excluded[i][j] ? Double.MAX_VALUE : cost[i][j] + node.pi[i] + node.pi[j];
    }
    int firstNeighbor;
    int secondNeighbor;
    // find the two cheapest edges from 0
    if (costWithPi[0][2] < costWithPi[0][1]) {
      firstNeighbor = 2;
      secondNeighbor = 1;
    } else {
      firstNeighbor = 1;
      secondNeighbor = 2;
    }
    for (int j = 3; j < n; j++) {
      if (costWithPi[0][j] < costWithPi[0][secondNeighbor]) {
        if (costWithPi[0][j] < costWithPi[0][firstNeighbor]) {
          secondNeighbor = firstNeighbor;
          firstNeighbor = j;
        } else {
          secondNeighbor = j;
        }
      }
    }
    addEdge(node, 0, firstNeighbor);
    Arrays.fill(node.parent, firstNeighbor);
    node.parent[firstNeighbor] = 0;
    // compute the minimum spanning tree on nodes 1..n-1
    double[] minCost = costWithPi[firstNeighbor].clone();
    for (int k = 2; k < n; k++) {
      int i;
      for (i = 1; i < n; i++) {
        if (node.degree[i] == 0) break;
      }
      for (int j = i + 1; j < n; j++) {
        if (node.degree[j] == 0 && minCost[j] < minCost[i]) i = j;
      }
      addEdge(node, node.parent[i], i);
      for (int j = 1; j < n; j++) {
        if (node.degree[j] == 0 && costWithPi[i][j] < minCost[j]) {
          minCost[j] = costWithPi[i][j];
          node.parent[j] = i;
        }
      }
    }
    addEdge(node, 0, secondNeighbor);
    node.parent[0] = secondNeighbor;
    node.lowerBound = Math.rint(node.lowerBound);
  }

  private void addEdge(Node node, int i, int j) {
    double q = node.lowerBound;
    node.lowerBound += costWithPi[i][j];
    node.degree[i]++;
    node.degree[j]++;
  }
}

class Node {
  public boolean[][] excluded;
  // Held--Karp solution
  public double[] pi;
  public double lowerBound;
  public int[] degree;
  public int[] parent;
}

class NodeComparator implements Comparator<Node> {
  public int compare(Node a, Node b) {
    return Double.compare(a.lowerBound, b.lowerBound);
  }
}
2020-07-28