Travelling Salesperson Problem (TSP)

• Wir haben ein Netzwerk von Städten, die wir besuchen wollen.
• Wir starten in einer Stadt und wollen alle anderen Städte besuchen.
• Wir brauchen jede Stadt nur einmal besuchen und müssen am Ende wieder in der Startstadt ankommen.
• Wir wollen die kürzeste Route finden.

Beispiel: Die Odyssee

Odysseus muss 22 Orte im Mittelmeerraum besuchen

      latlon = [
          (38.24, 20.42),
          (39.57, 26.15),
          (40.56, 25.32),
          (36.26, 23.12),
          (33.48, 10.54),
          (37.56, 12.19),
          (38.42, 13.11),
          (37.52, 20.44),
          (41.23, 09.10),
          (41.17, 13.05),
          (36.08, -5.21),
          (38.47, 15.13),
          (38.15, 15.35),
          (37.51, 15.17),
          (35.49, 14.32),
          (39.36, 19.56),
          (38.09, 24.36),
          (36.09, 23.00),
          (40.44, 13.57),
          (40.33, 14.15),
          (40.37, 14.23),
          (37.57, 22.56),
      ]
        

Modellieren können wir das als Graph $G=(V, E)$.

• Jeder Knoten $i \in V$ ist ein Ort.
• Die Kantenkosten $c_{(i, j)}$ sind die Distanzen.

• Im Odysseus Beispiel können wir als Distanz die Luftlinie zwischen den Orten nehmen.
• Hier wäre der dann Graph vollständig verbunden.

          from geopy.distance import distance

          distance(latlon[0], latlon[1]).km
        

Welche Variablen machen Sinn, um das als IP zu modellieren?

Wie bei anderen Graphenproblemen eine binäre Variable je Kante: \[ x_{ij} \in \{0, 1\} \quad \forall (i, j) \in E \]

Für jede Kante kennen wir die Distanz, wir wollen also die Gesamtdistanz minimieren.

\[ \min \sum_{(i, j) \in E} c_{ij} x_{ij} \]

Welche Constraints brauchen wir?

Jeder Knoten wird genau einmal verlassen:

\[ \sum_{(i,j) \in E} x_{ij} = 1 \quad \forall i \in V \]

Jeder Knoten wird genau einmal betreten:

\[ \sum_{(j, i) \in E} x_{ji} = 1 \quad \forall i \in V \]

Sind wir damit schon fertig?

          import mip, itertools

          model = mip.Model()

          model.objective.sense = mip.MINIMIZE
          variables = {}
          for p1, p2 in itertools.product(latlon, repeat=2):
              if p1 != p2:
                  variables[p1, p2] = model.add_var(var_type=mip.BINARY, obj= distance(p1, p2).km)

          for p in latlon:
              outflow = sum(variables[p, p2] for p2 in latlon if p != p2)
              model.add_constr(outflow == 1)
              inflow = sum(variables[p2, p] for p2 in latlon if p != p2)
              model.add_constr(inflow == 1)

          model.optimize()
        

Die Lösung besteht aus kleinen Subtouren.

Für eine korrekt Lösung müssen wir also noch etwas machen.

Idee:Zähle für jeden Knoten, an welcher Stelle der Route er kommt.

• Wir wählen einen beliebigen Knoten $v$ als den ersten aus.
• Jeder andere Knoten bekommt eine Variable, welche angibt an welcher Stelle der Knoten kommt:
  \[ y_i \geq 0 \quad \forall i \in V \setminus \{v\} \]

Damit der Nachfolger von $v$ den richtigen Wert zugewiesen bekommt:

\[ y_i \geq x_{vi} \quad \forall i \in V \setminus \{v\} \]

Wenn $j$ der Nachfolger von $i$ ist, soll gelten:

Wir müssen diese Ungleichung aber ausschalten, falls $j$ gar nicht der Nachfolger von $i$ ist:

Keine Variable kann einen Wert über $|V|$ bekommen, daher können wir $M = |V|$ setzen:

\[ y_j \geq y_i + 1 - |V| \cdot (1 - x_{ij}) \\[2mm] \forall i, j \in V \setminus \{v\}, i \neq j \]

\[ \begin{align*} \min \quad & \sum_{(i, j) \in E} c_{ij} x_{ij} \\ \text{s.t} \quad & \sum_{(i,j) \in E} x_{ij} = 1 && \forall i \in V \\ & \sum_{(j, i) \in E} x_{ji} = 1 && \forall i \in V \\ & y_j \geq y_i + 1 - |V| \cdot (1 - x_{ij}) && \forall i, j \in V \setminus \{v\}, i \neq j \\ & x_{ij} \in \{0, 1\} && \forall (i, j) \in E \\ & y_i \geq 0 && \forall i \in V \setminus \{v\} \end{align*} \]

          order_vars = {p: model.add_var(lb=0) for p in latlon[1:]}
          for p1, p2 in itertools.product(latlon, latlon[1:]):
              if p1 == latlon[0]:
                  model.add_constr( order_vars[p2] - variables[p1,p2] >= 0 )
              elif p1 != p2:
                  model.add_constr( order_vars[p2] >= order_vars[p1] + 1 - len(latlon) * (1 - variables[p1,p2]) )
        

Jetzt haben wir die korrekte Lösung.

Nachteil dieses Modells: Die Big-M Constraints machen das Modell langsam.

Finden wir bessere Nebenbedingungen, um Subtouren auszuschließen?

Wieviele ausgewählte Kanten gibt es innerhalb der Knoten einer Subtour?

Wieviele hätte es geben dürfen?

Für jede Knotenteilmenge muss gelten:

\[ \sum_{\substack{(i, j) \in E\\ i \in S, j \in S}} x_{ij} \leq |S| - 1 \quad \forall S \subset V, |S| \geq 2 \]

Das sind gute Nebenbedingungen.

Aber leider zu viele:
Anzahl an Teilmengen von $V$ ist $2^{|V|}$

Wir können schonmal versuchen alle 2er Subtouren auszuschließen:

\[ x_{i,j} + x_{j,i} \leq 1 \quad \forall i, j \in V, i \neq j \]

      for p1, p2 in itertools.product(latlon, repeat=2):
          if p1 != p2:
              model.add_constr(variables[p1, p2] + variables[p2, p1] <= 1)
        

Die Lösung ist zumindest besser als vorher mit 2er Subtouren.

Wie verhindern wir größere Subtouren?

Idee: Generiere Contraints nur bei Bedarf.

Wir können dem Solver einen Callback mitgeben, der für eine aktuelle Lösung neue Constraints erzeugt.

Wir müssen uns überlegen, wie wir eine Subtour finden.

Welche Knoten sind die einer Subtour?

Die einer Seite eines Cuts mit Wert kleiner 1.

Wir suchen also einen minimalen Schnitt (egal zwischen welchen Knoten) mit Wert kleiner 1.

Eine Strategie: Suche minimalen s-t-Cut zwischen
$v$ und allen anderen Knoten.

Vorteil: den Schnitt können wir sogar finden, wenn die Lösung noch fraktional ist.

Der ConstrsGenerator defininiert den Callback:

      class SubTourCutGenerator(mip.ConstrsGenerator):
          def generate_constrs(self, node_m: mip.Model, depth: int = 0, npass: int = 0):
              t_vars = node_m.translate(variables)
        

Durch Preprocessing werden die Variablen im Knoten geändert und müssen zurück übersetzt werden.

Finden einer neuen Constraint:

          g = nx.DiGraph()
          for (p1, p2), v in t_vars.items():
              g.add_edge(p1, p2, capacity=v.x)

          s = latlon[0]
          for t in latlon[1:]:
              val, (S, NS) = nx.minimum_cut(g, s, t)

              if val < 0.99:
                  node_m.add_constr(
                    sum(t_vars[p1, p2] for p1, p2 in itertools.product(S, repeat=2) if p1 != p2) <= len(S) - 1)
                  node_m.add_constr(
                    sum(t_vars[p1, p2] for p1, p2 in itertools.product(NS, repeat=2) if p1 != p2) <= len(NS) - 1)
        

Jetzt müssen wir den Callback noch im Modell registrieren.

          model.cuts_generator = SubTourCutGenerator()
          model.lazy_constrs_generator = SubTourCutGenerator()
        

Auch mit dem ConstrsGenerator macht es Sinn, die 2er Subtouren direkt auszuschließen

Jeder Callback Aufruf kostet Zeit, so sparen wir uns einige davon.

Jetzt erhalten wir eine optimale Lösung und sind deutlich schneller als mit den Hilfsvariablen.

Es gibt zwei Typen von erzeugten Constraints.

• Cuts: cuts_generator
• Lazy Constraints: lazy_constrs_generator

Cuts

Optionale Constraints, welche die LP Relaxierung besser machen.

Wird mit fraktionalen Lösungen aufgerufen.

Lazy Constraints

Relevante Constraints, ohne die eine Lösung unzulässig sein könnte.

Wir mit ganzzahligen Lösungen aufgerufen.

Nachteil von Lazy Constraints:

Der Solver kann weniger Preprocessing machen
Er kennt ja nicht das ganze Modell.

Kostet potentiell Performance.

In unserem Fall können wir das beste aus beiden Welten haben.

• Hilfsvariablen, damit das Modell von Anfang an vollständig ist.
• Cut Generator, um bei Bedarf bessere Nebenbedingungen zu erzeugen.
• Keine Lazy Constraints.

Mit diesem Code können wir schon recht große TSPs in akzeptabler Zeit optimal lösen:

Constraints bei Bedarf zu generieren kann oft Performance bringen.

Auch wenn die Anzahl Nebenbedingungen nicht exponentiell sondern z.B. nur kubisch ist.

Es gibt ein analoges Konzept, um Variablen nur bei Bedarf zu generieren:

Column Generation

Das ist aber ein gutes Stück umständlicher.