Оптимизированные алгоритмы TSP

Ответ 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);
  }
}

Ответ 2

Если ваш график удовлетворяет неравенству треугольника и вам нужна гарантия 3/2 в пределах оптимального, я предлагаю алгоритм christofides. Я написал реализацию в php на phpclasses.org.

Ответ 3

Начиная с 2013 года, можно решить для 100 городов, используя только точную формулировку в Cplex. Добавьте уравнения степени для каждой вершины, но включите ограничители, ограничивающие только то, что они появляются. Большинство из них не являются необходимыми. На этом примере Cplex.

Вы должны быть в состоянии решить 100 городов. Вам придется перебирать каждый раз, когда будет найден новый подтекст. Я привел пример здесь, а через пару минут и 100 итераций позже получил свои результаты.

Ответ 4

Я взял алгоритм Held-Karp из библиотеки Concorde, и 25 городов были решены за 0,15 секунды. Эта работа отлично подходит для меня! Вы можете извлечь код (написанный в ANSI C) из удерживаемого карпа из библиотеки concorde: http://www.math.uwaterloo.ca/tsp/concorde/downloads/downloads.htm. Если загрузка имеет расширение gz, это должно быть tgz. Возможно, вам придется переименовать его. Затем вы должны сделать небольшую корректировку для входа в VС++. Сначала возьмите файл holdkarp h и c (переименуйте его cpp) и другие около 5 файлов, внесите корректировки и он должен работать с вызовом CCheldkarp_small (...) с edgelen: euclid_ceiling_edgelen.

Ответ 5

TSP - проблема NP-hard. (Насколько нам известно) нет алгоритма NP-жестких задач, который выполняется в полиномиальное время, поэтому вы запрашиваете то, чего не существует.

Это либо достаточно быстро, чтобы закончить в разумные сроки, а затем это не точно или точно, но не закончится в вашей жизни за 100 городов.

Ответ 6

Чтобы дать глупый ответ: я тоже. В таком алгоритме все заинтригованы, но, как уже отмечали другие, я пока не существует (пока?). Esp ваша комбинация точных, 200 узлов, нескольких секунд выполнения и всего 200 строк кода невозможна. Вы уже знаете, что это NP трудно, и если вы получили небольшое представление об асимптотическом поведении, вы должны знать, что нет никакого способа достичь этого (за исключением того, что вы докажете, что NP = P, и даже я бы сказал, что это невозможно). Даже точные коммерческие решатели нуждаются в таких случаях гораздо больше, чем несколько секунд, и, как вы можете себе представить, у них гораздо больше 200 строк кода (даже если вы просто рассматриваете их ядра).

EDIT: Алгоритмы wiki являются "обычными подозреваемыми" в поле: Линейное программирование и ответвление. Их решения для экземпляров с тысячами узлов потребовали Years для решения (они просто сделали это с очень большим количеством процессоров параллельно, чтобы они могли делать это быстрее). Некоторые даже используют для специфических знаний о ветвях и границах для ограничения, поэтому они не являются общими подходами.

В ветке и связанной строке перечислены все возможные пути (например, с обратным отслеживанием) и применяется, когда у него есть решение, чтобы остановить начальную рекурсию, когда она может доказать, что результат не лучше, чем уже найденное решение (например, если вы просто посетили 2 из ваших городов, и путь уже длиннее, чем 200 экскурсия по городу. Вы можете отказаться от всех туров, которые начинаются с этой 2-х городской комбинации). Здесь вы можете вкладывать очень много особых знаний в функцию, которая говорит вам, что путь не будет лучше, чем уже найденное решение. Чем лучше, тем меньше путей вы должны смотреть, тем быстрее ваш алгоритм.

Линейное программирование - это метод оптимизации, поэтому решайте задачи линейного неравенства. Он работает в полиномиальное время (симплекс только практически, но это не имеет значения здесь), но решение является реальным. Когда у вас есть дополнительное ограничение, что решение должно быть целым, оно получает NP-complete. В небольших случаях это возможно, например. один метод его решения, затем посмотрите, какая переменная решения нарушает целочисленную часть и добавляет дополнительные неравенства для ее изменения (это называется плоскостью резания, имя указывает на то, что неравенства определяют (более объемную) плоскость, пространство решений является многогранником и, добавляя дополнительные неравенства, вы что-то вырезаете плоскостью из политопа). Тема очень сложная и даже простой простой симплекс трудно понять, когда вы не хотите глубоко погрузиться в математику. Есть несколько хороших книг, один из лучших - от Chvatal, Linear Programming, но есть еще несколько.

Ответ 7

У меня есть теория, но у меня никогда не было времени ее преследовать:

TSP - это ограниченная проблема (одна форма, где все точки лежат по периметру), где оптимальное решение - это решение, имеющее самый короткий периметр.

Существует множество простых способов получить все точки, которые лежат на минимальном граничном периметре (представьте себе, что большая эластичная лента растянута вокруг кучки гвоздей на большой доске).

Моя теория заключается в том, что если вы начнете надавливать на эластичную ленту так, чтобы длина полосы увеличивалась на ту же величину между соседними точками по периметру, и каждый сегмент оставался в форме эллиптической дуги, растянутый эластичный пересекают точки на оптимальном пути до точек пересечения на неоптимальных дорожках. См. эту страницу на mathopenref.com на эллипсах рисования - особенно шаги 5 и 6. Точки на ограничивающем периметре можно рассматривать как координаты эллипс (F1, F2) на изображениях ниже.

То, что я не знаю, - это если процесс "растяжения пузырьков" должен быть reset после добавления каждой новой точки или если существующие "пузыри" продолжают расти, и каждая новая точка по периметру вызывает только локализованный "пузырь", чтобы превратиться в два отрезка. Я оставлю это для вас, чтобы понять.