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 valuelet bestsolsofar = nulllet cursol be the partial solution with no variables assignedloop 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 cursolend loop分支定界的想法是,存在一个部分解决方案的搜索树。解决Held–Karp的要点是LP的值最多为最佳行程的OPT长度,但推测至少为3/4
OPT(实际上,通常更接近OPT)。
我遗漏的伪代码中的一个细节是如何选择分支变量。通常,目标是首先做出“艰难”的决定,因此,修复其值已经接近0或1的变量可能是不明智的。一种选择是选择最接近0.5的值,但还有很多其他的选择。
编辑
Java实现。198条非空白,非注释行。我忘记了1树不能将变量分配给1,所以我通过查找1树的度大于2的顶点进行分支并依次删除每个边。该程序接受TSPLIB实例的
EUC_2D形式,例如,
eil51.tsp与
eil76.tsp和
eil101.tsp和
lin105.tsp从http://www2.iwr.uni-
heidelberg.de/groups/comopt/software/TSPLIB95/tsp/。
// simple exact TSP solver based on branch-and-bound/Held--Karpimport 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("%ft%ft%ft%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); }}


