Ответ 1
200 строк и никакие библиотеки не являются жестким ограничением. Передовые решатели используют ветку и связаны с релаксацией Хелд-Карпа, и я не уверен, что даже самая базовая версия этой системы вписывается в 200 нормальных линий. Тем не менее, здесь схема.
Сопряженный Карп
Один из способов записи TSP в виде целочисленной программы выглядит следующим образом (Данциг, Фулкерсон, Джонсон). Для всех ребер e константа w e обозначает длину ребра e, а переменная x e равна 1, если ребро e находится в туре и 0 в противном случае. Для всех подмножеств S вершин, ∂ (S) обозначает ребра, соединяющие вершину в S с вершиной, не принадлежащей S.
минимизировать сумму ребра e w e x e
с учетом
1. для всех вершин v, sum ребра e в ∂ ({v}) x e= 2
2. для всех непустых собственных подмножеств S вершин sum ребра e в ∂ (S) x e ≥ 2
3. для всех ребер e в E, x e в {0, 1}
Условие 1 гарантирует, что множество ребер представляет собой набор туров. Условие 2 гарантирует наличие только одного. (В противном случае пусть S - множество вершин, посещенных одним из туров.) Релаксация Held-Karp получается путем внесения этого изменения.
3. для всех ребер e в E, x e в {0, 1}
3. для всех ребер e в E, 0 ≤ x e ≤ 1
Held-Karp является линейной программой, но имеет экспоненциальное число ограничений. Один из способов его решения - ввести множители Лагранжа, а затем сделать субградиентную оптимизацию. Это сводится к циклу, который вычисляет минимальное связующее дерево, а затем обновляет некоторые векторы, но детали могут быть задействованы. Помимо "Held-Karp" и "subgradient (спуск | оптимизация)", "1-дерево" - еще один полезный термин поиска.
(Более медленная альтернатива заключается в том, чтобы записывать решатель 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 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.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--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);
}
}
