algorithm - travelling - Algoritmos TSP optimizados

problema del agente viajero tsp travelling salesman problem

Question

Me interesan las formas de mejorar o encontrar algoritmos que puedan resolver el problema del vendedor itinerante por aproximadamente n = 100 to 200 ciudades.

El enlace de wikipedia que presenté enumera varias optimizaciones, pero lo hace a un nivel bastante alto, y no sé cómo implementarlas en el código.

Hay solucionadores de resistencia industrial, como Concorde , pero son demasiado complejos para lo que quiero, y las soluciones clásicas que inundan las búsquedas de TSP presentan algoritmos aleatorios o algoritmos clásicos de retroceso o programación dinámica que solo funcionan para 20 ciudades.

Entonces, ¿alguien sabe cómo implementar un simple (simplemente quiero decir que una implementación no toma más de 100-200 líneas de código) solucionador de TSP que funciona en un tiempo razonable (unos pocos segundos) para al menos 100 ciudades? Solo estoy interesado en soluciones exactas.

Puede suponer que la entrada se generará de forma aleatoria, por lo que no me importan las entradas destinadas específicamente a romper un determinado algoritmo.

Answer 1

200 líneas y ninguna biblioteca es una restricción difícil. Los solucionadores avanzados usan branch y bound con la relajación de Held-Karp, y no estoy seguro de que incluso la versión más básica de eso encaje en 200 líneas normales. Sin embargo, aquí hay un esquema.

Held Karp

Una forma de escribir TSP como un programa entero es la siguiente (Dantzig, Fulkerson, Johnson). Para todos los bordes e, la constante w _e denota la longitud del borde e, y la variable x _e es 1 si el borde e está en el recorrido y 0 en caso contrario. Para todos los subconjuntos S de vértices, ∂ (S) denota los bordes que conectan un vértice en S con un vértice no en S.

minimizar los _{bordes de} suma _e w _e x _e
sujeto a
1. para todos los vértices v, suma los _{bordes e en ∂ ({v})} x _e = 2
2. para todos los subconjuntos propios no vacíos S de vértices, _bordes suma _{e en ∂ (S)} x _e ≥ 2
3. para todos los bordes e en E, x _e en {0, 1}

La condición 1 asegura que el conjunto de bordes sea una colección de recorridos. La condición 2 asegura que solo hay una. (De lo contrario, que S sea el conjunto de vértices visitados por uno de los recorridos). La relajación de Held-Karp se obtiene haciendo este cambio.

3. para todos los bordes e en E, x _e en {0, 1}
3. para todos los bordes e en E, 0 ≤ x _e ≤ 1

Held-Karp es un programa lineal pero tiene un número exponencial de restricciones. Una forma de resolverlo es introducir multiplicadores de Lagrange y luego realizar una optimización subgradiente. Eso se reduce a un bucle que calcula un árbol de expansión mínimo y luego actualiza algunos vectores, pero los detalles están más o menos involucrados. Además de "Held-Karp" y "subgradiente (descenso | optimización)", "1-árbol" es otro término de búsqueda útil.

(Una alternativa más lenta es escribir un solucionador de LP e introducir restricciones de subterfugio, ya que se violan por optima anterior. Esto significa escribir un solucionador de LP y un procedimiento de corte mínimo, que también es más código, pero podría extenderse a un TSP más exótico restricciones)

Sucursal y atado

Por "solución parcial", me refiero a una asignación parcial de variables a 0 o 1, donde un borde asignado 1 está definitivamente en el recorrido, y un borde asignado 0 definitivamente está fuera. La evaluación de Held-Karp con estas restricciones laterales brinda un límite inferior en el recorrido óptimo que respeta las decisiones ya tomadas (una extensión).

Branch and bound mantiene un conjunto de soluciones parciales, al menos una de las cuales se extiende a una solución óptima. El pseudocódigo para una variante, la búsqueda en profundidad con el mejor retroceso inicial es el siguiente.

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

La idea de branch y bound es que hay un árbol de búsqueda de soluciones parciales. El objetivo de resolver Held-Karp es que el valor del LP es como máximo la longitud OPT del recorrido óptimo, pero también se supone que es al menos 3/4 OPT (en la práctica, generalmente más cerca de OPT).

El único detalle en el pseudocódigo que he omitido es cómo elegir la variable de bifurcación. El objetivo suele ser tomar las decisiones "difíciles" primero, por lo que es probable que corregir una variable cuyo valor ya está cerca de 0 o 1 no sea prudente. Una opción es elegir la más cercana a 0.5, pero hay muchas, muchas otras.

EDITAR

Implementación de Java. 198 líneas no en blanco, sin comentario. Olvidé que los árboles 1 no funcionan con la asignación de variables a 1, así que me ramifico al encontrar un vértice cuyo 1-árbol tiene un grado> 2 y elimino cada borde por turno. Este programa acepta instancias TSPLIB en formato EUC_2D , por ejemplo, eil51.tsp y eil76.tsp y eil101.tsp y lin105.tsp desde 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);
  }
}

Answer 2

A partir de 2013, es posible resolver para 100 ciudades utilizando solo la formulación exacta en Cplex. Agregue ecuaciones de grados para cada vértice, pero incluya restricciones para evitar las sutilezas solo cuando aparezcan. La mayoría de ellos no son necesarios. Cplex tiene un ejemplo sobre esto.

Deberías poder resolver por 100 ciudades. Tendrá que repetir cada vez que encuentre una nueva suburbanidad. Corrí un ejemplo aquí y en un par de minutos y 100 iteraciones más tarde obtuve mis resultados.

Answer 3

Para dar una respuesta tonta: yo también. Todo el mundo está interesado en dicho algoritmo, pero como otros ya han dicho: no existe (¿todavía?). Esp su combinación de 200 nodos exactos, tiempo de ejecución de pocos segundos y solo 200 líneas de código son imposibles. Ya sabes que es NP difícil y si tienes la menor impresión de comportamiento asintótico debes saber que no hay forma de lograrlo (excepto que demuestres que NP = P, e incluso eso diría que eso no es posible). Incluso los solucionadores comerciales exactos necesitan esos casos mucho más que algunos segundos y, como se puede imaginar, tienen más de 200 líneas de código (incluso cuando solo considera sus kernels).

EDITAR: Los algoritmos de la wiki son los "sospechosos habituales" del campo: Programación lineal y sucursal. Sus soluciones para las instancias con miles de nodos tardaron años en resolverse (simplemente lo hicieron con muchas CPUs paralelas, por lo que pueden hacerlo más rápido). Algunos incluso usan para el conocimiento específico del problema de ramificación y delimitación para el límite, por lo que no son enfoques generales.

Branch and bound solo enumera todas las rutas posibles (por ejemplo, con backtracking) y aplica una vez que tiene una solución para detener una recursión iniciada cuando puede demostrar que el resultado no es mejor que la solución ya encontrada (por ejemplo, si acaba de visitar 2 de sus ciudades y el camino ya es más largo que 200 recorridos por la ciudad. Puede descartar todos los tours que comiencen con esa combinación de 2 ciudades). Aquí puede invertir mucho conocimiento específico del problema en la función que le indica que la ruta no va a ser mejor que la solución ya encontrada. Cuanto mejor es, menos caminos hay que mirar, más rápido es su algoritmo.

La programación lineal es un método de optimización, así que resuelve problemas de desigualdad lineal. Funciona en tiempo polinomial (simplex solo de forma práctica, pero eso no importa aquí), pero la solución es real. Cuando tiene la restricción adicional de que la solución debe ser entera, obtiene NP completa. Para pequeñas instancias es posible, por ejemplo, un método para resolverlo, luego buscar qué variable de la solución viola la parte entera y agregar desigualdades de suma para cambiarlo (esto se llama plano de corte, el nombre proviene del hecho de que las desigualdades definen (plano de mayor dimensión), el espacio de la solución es un polytop y al agregar desigualdades adicionales, se corta algo con un plano del polytop). El tema es muy complejo e incluso un símplex simple general es difícil de entender cuando no quieres profundizar en las matemáticas. Hay varios buenos libros sobre, uno de los mejores es de Chvatal, Programación lineal, pero hay varios más.

Answer 4

Si su gráfica satisface la desigualdad del triángulo y desea una garantía de 3/2 dentro del óptimo, sugiero el algoritmo de christofides. Escribí una implementación en php en phpclasses.org.

Answer 5

TSP es un problema NP-difícil. (Por lo que sabemos) no existe un algoritmo para problemas NP-hard que se ejecuta en tiempo polinomial, por lo que pide algo que no existe.

Es lo suficientemente rápido como para terminar en un tiempo razonable y luego no es exacto o exacto, pero no terminará en tu vida en 100 ciudades.

Answer 6

Tengo una teoría, pero nunca tuve tiempo de buscarla:

El TSP es un problema límite (forma única donde todos los puntos se encuentran en el perímetro) donde la solución óptima es aquella que tiene el perímetro más corto.

Hay muchas formas simples de obtener todos los puntos que se encuentran en un perímetro límite mínimo (imagina una gran banda elástica estirada alrededor de un montón de uñas en una tabla grande).

Mi teoría es que si comienzas a empujar la banda elástica para que la longitud de la banda aumente en la misma cantidad entre puntos adyacentes en el perímetro, y cada segmento permanezca en forma de arco elíptico, el elástico estirado cruzará puntos en el camino óptimo antes de cruzar puntos en caminos no óptimos. Consulte esta página en mathopenref.com para dibujar elipses, especialmente los pasos 5 y 6. Los puntos en el perímetro delimitador se pueden ver como puntos focales de la elipse (F1, F2) en las imágenes a continuación.

Lo que no sé es si el proceso de "estiramiento de burbujas" debe restablecerse después de agregar cada nuevo punto, o si las "burbujas" existentes continúan creciendo y cada nuevo punto en el perímetro provoca solo la "burbuja" localizada. convertir en dos segmentos de línea. Lo dejo para que lo descubras.

Answer 7

Tomé el algoritmo Held-Karp de la biblioteca concorde y 25 ciudades se resuelven en 0.15 segundos. ¡Este rendimiento es perfectamente bueno para mí! Puede extraer el código (escrito en ANSI C) de held-karp de la biblioteca concorde: http://www.math.uwaterloo.ca/tsp/concorde/downloads/downloads.htm . Si la descarga tiene la extensión gz, debería ser tgz. Es posible que necesite cambiarle el nombre. Luego debe hacer pequeños ajustes para ingresar en VC ++. Primero tome el archivo holdkarp h y c (cámbiele el nombre cpp) y otros aproximadamente 5 archivos, haga los ajustes necesarios y llame al CCheldkarp_small (...) con edgelen: euclid_ceiling_edgelen.