Problem komiwojażera w praktyce

Korzystając z popularnych serwisów internetowych, błyskawicznie znajdziemy najkrótszą trasę między dwoma miastami. A co, gdybyśmy chcieli znaleźć najkrótszą trasę, która pozwoli po wyruszeniu z domu odwiedzić wszystkie interesujące nas miasta i wrócić do punktu wyjścia?

Możemy to pytanie sformalizować w następujący sposób. Wybrane miasta (wraz z tym, w którym mieszkamy) numerujemy od $|1$ do $|n.$ Zakładamy, że mamy tabelę $n × n$ o wartościach $|di j$ dla $i, j = 1,...,n,$ gdzie $|d i j$ oznacza odległość (np. podaną przez nasz ulubiony serwis) z miasta $|i$ do miasta $j.$ Tabela jest symetryczna, tzn. $|di j = d ji.$ Naszym celem jest znalezienie kolejności odwiedzania miast, czyli permutacji $|v1,...,vn,$ która minimalizuje długość trasy, tzn. $|Pni 0dvivi 1,$ przyjmując $v0 = vn.$ To zadanie znane jest jako problem komiwojażera i od kilkudziesięciu lat spędza sen z oczu informatykom.

Z jednej strony jest to problem NP-trudny, a najlepszy w sensie złożoności pesymistycznej algorytm działa w czasie $n 2 O(2 n ).$ Z drugiej strony, istnieją implementacje, które dla instancji pochodzących ze świata rzeczywistego uzyskują spektakularne wyniki (tabela poniżej), ale ich pesymistyczna złożoność czasowa nie jest znana.

Wszystkie powyższe rekordy dotyczą optymalnych rozwiązań dla rzeczywistych map połączeń drogowych (w USA, Niemczech, Szwecji), z wyjątkiem rekordu Manfreda Padberga i Giovanniego Rinaldiego, którzy znaleźli najkrótszą trasę dla wiertarki wiercącej otwory w płytce obwodu drukowanego. Aktualnie rozwiązanie problemu na zwykłym laptopie dla 1000 miast to kwestia kilkudziesięciu sekund, ale już wynik z ostatniego wiersza tabelki (2016) uzyskano za pomocą 310 procesorów pracujących przez 8 miesięcy. Choć pierwszy wiersz tabelki wygląda skromnie, w rzeczywistości stanowi wielkie osiągnięcie, a kolejne rekordy były wynikiem rozwinięcia technik wypracowanych przez George'a Dantziga, Delberta Fulkersona i Selmera Johnsona. W dalszej części artykułu zobaczymy, jak działa ich metoda na nieco bliższym nam przykładzie 51 dużych miast w Polsce.

Kluczowy pomysł polega na zapisaniu naszego problemu jako zbioru równań i nierówności liniowych. Dla każdej pary różnych miast $|i, j = 1,...,51$ wprowadźmy zmienną $|x , i j$ przy czym $x i j$ oraz $x ji$ traktujemy jako dwie nazwy tej samej zmiennej. Na początek przyjmijmy, że zmienne $|xi j$ mogą przyjmować tylko dwie wartości: 0 i 1. Myślimy o nich w ten sposób, że połączenie między $i$ oraz $j$ jest używane w naszym rozwiązaniu wtedy i tylko wtedy, gdy $xi j = 1.$ Zauważmy, że długość takiej trasy wynosi $|P1Di@ jDndi jxi j.$

Czy możemy za pomocą prostych równań i nierówności wymusić, aby spełniające je wartości zmiennych odpowiadały dokładnie możliwym trasom komiwojażera (czyli cyklom długości $|n$ )? Pierwsza obserwacja jest prosta: dla każdego miasta $i$ dokładnie dwie zmienne $|xi j$ powinny mieć wartość 1, co można wyrazić jako $|Pj xixi j = 2.$ To jednak nie wystarczy, gdyż zamiast jednej trasy moglibyśmy dostać kilka krótszych, rozłącznych cykli (w skrajnym przypadku $|n/3$ rozłącznych trójkątów).

Aby zagwarantować, że miasta z pewnego zbioru $S$ nie tworzą krótszej trasy, możemy dodać warunek eliminacji podtras, który mówi, że do tego zbioru musimy co najmniej jeden raz wejść i co najmniej jeden raz z niego wyjść. Zwykle ten warunek zapisujemy w równoważnej formie $|P x ⩽ S − 1 i, j>S,i@ j i j$ (czyli $P x ⩾ 2 i>S, j~>S i j$ ). Równoważność łatwo dostaniemy, mnożąc oryginalną nierówność przez $|− 1,$ dodając do niej równości $|Pj xixi j = 2$ dla wszystkich $i∈ S$ i dzieląc otrzymaną nierówność przez $|2.$ Niestety, taki warunek musimy dodać dla każdego zbioru miast $|S$ o mocy między $|3$ a $⌊n/2⌋.$ Otrzymujemy następujące zadanie, które jest przykładem tzw. programu liniowego całkowitoliczbowego.

Zminimalizuj $| Q di jxi j 1Di@ jDn$ przy ograniczeniach

$|P jxixi j = 2,$ dla $i = 1,...,n,$
$|Pi, j> S,i@ jxi j⩽ S − 1,$ dla $|S ⊆{1,...,n}$ i $n 3 ⩽ S ⩽⌊ 2⌋ ,$
$xi j∈ {0,1},$ dla $1⩽ i < j ⩽n.$

Skąd pomysł na zapisanie problemu w języku równań i nierówności liniowych? Otóż teraz rozluźniamy założenia naszego zadania: warunki $|xi j ∈{0,1}$ zamieniamy na $|0⩽ xi j⩽ 1$ i szukamy rozwiązania w liczbach rzeczywistych. Wówczas otrzymujemy tzw. program liniowy. Tak się składa, że nasi bohaterowie umieli szybko rozwiązywać programy liniowe, gdyż jeden z nich, George Dantzig, w 1947 roku opublikował do dziś najskuteczniejszy w praktyce algorytm realizujący ten cel, nazywany algorytmem simplex. Jest on obecnie dostępny w bibliotekach dla wielu popularnych języków programowania, np. przygotowując ten artykuł, użyliśmy modułu PuLP dla języka Python.

Plan jest następujący. Jeśli szczęśliwie jako rozwiązanie optymalne naszego programu liniowego dostaniemy same zera i jedynki, to będą one definiowały pewną trasę, która musi być optymalna, ponieważ każda trasa komiwojażera spełnia wszystkie warunki programu. Niestety może się okazać, że dostaniemy rozwiązanie ułamkowe, w którym niektóre zmienne będą miały wartość w przedziale $(0,1).$ Wtedy będziemy musieli jakoś zareagować.

Dodatkowy kłopot polega na tym, że nasz program liniowy ma gigantyczną liczbę warunków eliminacji podtras: tylko dla $S = 25$ jest ich $|(5125) = 247959266474052.$ To za dużo nawet dla współczesnych komputerów. Pomysł Dantziga i jego kolegów polegał na tym, aby zacząć od prostego programu liniowego postaci:

Zminimalizuj $| Q d x 1Di@ jDn i j i j$ przy ograniczeniach

$|P jxixi j = 2,$ dla $|i = 1,...,n,$
$0 ⩽ x ⩽ 1, i j$ dla $1⩽ i < j ⩽ n$ ;

a następnie dodawać warunki w miarę potrzeby. Taki program ma $|(51) = 1275 2$ zmiennych i 2601 warunków liniowych. Uruchamiamy na komputerze algorytm simplex i w ułamku sekundy dostajemy poniższe rozwiązanie.

Liniami ciągłymi oznaczyliśmy zmienne o wartości $1,$ natomiast przerywanymi o wartości $|12$ (innych ułamkowych wartości w uzyskanym rozwiązaniu nie było). Widzimy, że możemy dodać warunki eliminacji podtras dla zbiorów ${14,33,41},{2,6,13},{32,34, 18},{17,24,42, 48}$ oraz $|{12,15,36}.$ Chociaż miasta $4,21$ i $49$ nie utworzyły krótkiego cyklu, zauważmy, że dla $S = {4,21,49}$ oraz dla $S = {11,44,19}$ warunek eliminacji podtras też nie jest spełniony; dodajemy więc jeszcze te dwa warunki i ponownie uruchamiamy algorytm.

W nowym rozwiązaniu mamy nowe krótkie cykle, dodajemy warunki eliminacji podtras dla kolejnych dwóch zbiorów $|{2,6,13,17,42, 24,48}$ oraz $|{4,21,0,35,45,37,33,14, 41,34, 18,32,49}$ i generujemy kolejne rozwiązanie.

Otrzymaliśmy dwa cykle, a więc dodajemy warunek eliminacji podtras dla krótszego z nich o zbiorze wierzchołków ${11,19,44,47}.$ Efekt kolejnego uruchomienia algorytmu simplex widzimy poniżej.

Niestety wszystkie warunki eliminacji podtras są spełnione, a jednak nie uzyskaliśmy pojedynczego cyklu. Jest to spowodowane dopuszczeniem ułamkowych wartości zmiennych. Ratunkiem w takiej sytuacji jest znalezienie nowej nierówności, którą spełniają trasy komiwojażera, ale nie spełnia jej nasze rozwiązanie. Takie nierówności nazywamy płaszczyznami odcinającymi. Dlaczego? Wyobraźmy sobie, że nasz program liniowy ma tylko dwie zmienne, a wszystkie warunki to nierówności. Wówczas zbiór punktów spełniających te warunki jest wielokątem na płaszczyźnie, jak poniżej.

Zminimalizuj $|x 1$ przy ograniczeniach

$x2 ⩽ x1 + 4,$
$2x +x ⩽ 8, 1 2$
$2x2 + x1⩽ 0,$
$x ⩽ 0. 2$

Algorytm simplex ma tę własność, że zwraca nam rozwiązanie optymalne, które jest wierzchołkiem wielokąta. Tymczasem my szukamy całkowitoliczbowego rozwiązania optymalnego. Możemy wówczas do programu liniowego dodać nierówność, która oddzieli nasz wierzchołek od wszystkich punktów o współrzędnych całkowitych wewnątrz wielokąta, odcinając fragment wielokąta za pomocą linii prostej (na rysunku fioletowa linia jest przykładem takiej prostej). W wyniku tej procedury możemy znów dostać rozwiązanie niecałkowitoliczbowe, wówczas powtarzamy operację. Z analogiczną sytuacją mamy do czynienia w większej liczbie wymiarów, gdzie uogólnieniem linii jest hiperpłaszczyzna.

Wróćmy do ostatniego rozwiązania i spójrzmy na trójkąt $|H = {19,20,50},$ którego krawędzie odpowiadają zmiennym połówkowym $|x1920,x2050$ i $x5019,$ oraz na sąsiadujące z nim połączenia $T = {1911,2031,5010}.$ W naszym rozwiązaniu suma $|s = x1920 + x2050 +x5019 + x1911 ++x2031 + x5010$ wynosi $1 4 2 .$ Tymczasem w dowolnej trasie komiwojażera $|s⩽ 4.$ Dlaczego? Z warunku $Pj xi xi j = 2$ dla $|i∈H$ wiemy, że

$pict$

Dodając stronami te nierówności i trywialną $Pe> T xe ⩽ T ,$ otrzymujemy $|2s⩽ 9.$ Ale jeśli zmienne opisują trasę komiwojażera, to $s$ jest liczbą całkowitą, a wtedy otrzymana nierówność jest równoważna $2s ⩽ 8,$ czyli $|s⩽ 4.$ Zachęcamy w tym miejscu Pracowitego Czytelnika, aby wykazał, że ogólniej, dla dowolnego zbioru miast $H$ i nieparzystej liczby połączeń $|T$ wychodzących z $|H,Pi, j>H,i@ j xi j + Pi j>T xi j⩽ H + STS−2 1.$ Są to tzw. nierówności 2-skojarzeniowe, gdyż zostały odkryte w kontekście problemu skojarzeń w grafach. Do naszego programu dodajmy nierówność 2-skojarzeniową $s⩽ 4,$ która pozwoli się pozbyć połówkowego trójkąta $|{19, 20,50}.$ Okazuje się, że w rozwiązaniu optymalnym takiego programu liniowego wszystkie zmienne mają już wartość 0 lub 1, a odpowiadające im rozwiązanie tworzy pojedynczy cykl, a więc jest optymalną trasą komiwojażera. Wygląda ona następująco.

Dantzig z kolegami nie odkryli nierówności 2-skojarzeniowych, podobnie jak w naszym przykładzie nierówności eliminacji podtras nie były wystarczające, więc musieli użyć dwóch dziwacznych nierówności. W kolejnych latach odkryto metody automatycznego znajdowania płaszczyzn odcinających, jednak wciąż najbardziej skuteczne okazują się rodziny płaszczyzn (czyli po prostu nierówności) projektowane z myślą o konkretnym problemie, jak nasze nierówności 2-skojarzeniowe. Dla problemu komiwojażera odkryto wiele takich rodzin. Dla niektórych z nich jesteśmy w stanie szybko (w czasie wielomianowym) znajdować nierówność niespełnioną przez dane rozwiązanie, o ile taka istnieje (jest to możliwe dla nierówności eliminacji podtras i nierówności 2-skojarzeniowych). Dla innych stosujemy heurystyki, które wprawdzie działają szybko, ale mogą pozostawić jakąś niespełnioną nierówność.

Dociekliwy Czytelnik na pewno zadaje sobie teraz pytanie, co się dzieje, gdy program rozwiązujący problem komiwojażera nie miał tyle szczęścia co my z mapą Polski i w pewnym kroku generuje rozwiązanie, które nie jest całkowitoliczbowe, a mimo to wszystkie rodziny nierówności, które jest w stanie sprawdzić, są spełnione. Odpowiedź jest prosta: wówczas używa się brutalnej siły. Wybieramy dowolną zmienną $xi j$ o wartości w przedziale $|(0,1)$ i rekurencyjnie rozwiązujemy dwa podproblemy. W pierwszym dodajemy warunek $|xi j = 1,$ w drugim $xi j = 0.$ Zauważmy, że każda optymalna trasa komiwojażera spełnia warunki dokładnie jednego z tych podproblemów. W każdym z podproblemów możemy znowu szukać płaszczyzn odcinających i rozgałęziać się na kolejne podproblemy. W ten sposób znajdziemy wiele tras komiwojażera, z których najkrótsza będzie optymalna.

Liczba podproblemów może szybko wymknąć się spod kontroli, bo $|k$ poziomów rozgałęzień to $|2k$ podproblemów. Z pomocą przychodzą wówczas heurystyki znajdujące niekoniecznie optymalne trasy komiwojażera (obecnie najskuteczniejsza, heurystyka Lina-Kernighana-Helsgauna, dla mapy świata zawierającej $|1904711$ miast znajduje trasę dłuższą od optymalnej o co najwyżej 0,023 %). Powiedzmy, że heurystyka znalazła nam trasę komiwojażera o długości $u,$ i rozważmy podproblem, w którym wartość rozwiązania optymalnego programu liniowego wynosi $|l.$ Każda trasa, którą możemy znaleźć w tym podproblemie, będzie miała długość co najmniej $l,$ a więc jeśli $l⩾ u,$ to taki podproblem możemy zignorować, gdyż nie pozwoli on znaleźć lepszej trasy niż ta, którą już mamy. Dzięki temu możliwe jest znaczne ograniczenie szybkiego wzrostu liczby podproblemów.

Okazuje się, że dla danych pochodzących ze świata rzeczywistego i dostatecznie bogatych rodzin płaszczyzn oddzielających liczba generowanych podproblemów nie jest astronomiczna. Jest to jednak jedynie obserwacja empiryczna. W 2018 roku Stefan Hougardy i Xianghui Zhong wymyślili sposób generowania sztucznych, wyjątkowo trudnych instancji problemu komiwojażera. Na podstawie testów oszacowali, że Concorde, najlepszy obecnie program korzystający z naszkicowanej w tym artykule metody, uruchomiony na współczesnym komputerze potrzebuje $3⋅1022$ lat obliczeń dla wygenerowanej przez nich instancji 1000 miast. To pokazuje, że przynajmniej na razie NP-trudności problemu nie dało się oszukać. Na szczęście duże instancje problemu komiwojażera rozwiązane optymalnie lub niemal optymalnie pokazują, że świat wokół nas rzadko jest aż tak złośliwy jak generator Hougardy'ego i Zhonga.

Optymalna trasa dla 51 polskich miast zaprezentowana w tym artykule została obliczona na podstawie tabeli $1275$ długości najkrótszych połączeń drogowych między nimi, znalezionych za pomocą serwisu Google maps we wrześniu 2019 roku. Dane oraz prosty skrypt w języku Python, który znajduje optymalną trasę i tworzy rysunki do tego artykułu, są dostępne na github.com/lkowalik/tsp. Zachęcamy do eksperymentów!