Jaki jest minimalny koszt połączenia wszystkich wysp?

Jest ruszt o rozmiarach N x m . Niektóre komórki to wyspy oznaczone „0”, a inne to woda . Każda komórka wodna ma na sobie liczbę oznaczającą koszt mostu wykonanego w tej komórce. Musisz znaleźć minimalny koszt, za jaki wszystkie wyspy mogą być połączone. Komórka jest połączona z inną komórką, jeśli ma wspólną krawędź lub wierzchołek.

Jakiego algorytmu można użyć do rozwiązania tego problemu? Czego można użyć jako metody brutalnej siły, jeśli wartości N, M są bardzo małe, powiedzmy NxM <= 100?

Przykład : Na podanym obrazku zielone komórki oznaczają wyspy, niebieskie pola - wodę, a jasnoniebieskie komórki - komórki, na których ma być wykonany most. Zatem dla poniższego obrazu odpowiedź będzie wynosić 17 .

Początkowo myślałem o oznaczeniu wszystkich wysp jako węzłów i połączeniu każdej pary wysp najkrótszym mostem. Wtedy problem mógłby zostać zredukowany do Minimalnego drzewa rozpinającego, ale w tym podejściu przegapiłem przypadek, w którym krawędzie zachodzą na siebie. Na przykład na poniższym obrazku najkrótsza odległość między dowolnymi dwiema wyspami wynosi 7 (zaznaczona na żółto), więc przy użyciu minimalnych drzew rozpinających odpowiedź będzie wynosić 14 , ale odpowiedź powinna wynosić 11 (zaznaczona na jasnoniebiesko).

Aby podejść do tego problemu, użyłbym struktury programowania liczb całkowitych i zdefiniowałbym trzy zestawy zmiennych decyzyjnych:

• x_ij : Zmienna wskaźnikowa binarna określająca, czy budujemy most na wodzie (i, j).
• y_ijbcn : Wskaźnik binarny określający, czy położenie wody (i, j) jest n ^ tą lokalizacją łączącą wyspę b z wyspą c.
• l_bc : Zmienna wskaźnika binarnego określająca, czy wyspy b i c są bezpośrednio połączone (aka można chodzić tylko po polach mostu od b do c).

W przypadku kosztów budowy mostów c_ij , docelową wartością do zminimalizowania jest sum_ij c_ij * x_ij. Do modelu musimy dodać następujące ograniczenia:

• Musimy upewnić się, że zmienne y_ijbcn są prawidłowe. Zawsze możemy dotrzeć do placu wodnego tylko wtedy, gdy zbudujemy tam most, a więc y_ijbcn <= x_ijdla każdego miejsca wodnego (i, j). Ponadto y_ijbc1musi wynosić 0, jeśli (i, j) nie graniczy z wyspą b. Wreszcie, dla n> 1, y_ijbcnmożna użyć tylko wtedy, gdy w kroku n-1 użyto sąsiedniej lokalizacji wodnej. Definiowanie N(i, j)jako kwadraty wody sąsiednie (i, j), jest równoważne z y_ijbcn <= sum_{(l, m) in N(i, j)} y_lmbc(n-1).
• Musimy upewnić się, że zmienne l_bc są ustawione tylko wtedy, gdy b i c są połączone. Jeśli zdefiniujemy I(c)lokalizacje graniczące z wyspą c, można to osiągnąć za pomocąl_bc <= sum_{(i, j) in I(c), n} y_ijbcn .
• Musimy zapewnić bezpośrednie lub pośrednie połączenie wszystkich wysp. Można to osiągnąć w następujący sposób: dla każdego niepustego właściwego podzbioru S wysp należy wymagać, aby co najmniej jedna wyspa w S była połączona z co najmniej jedną wyspą w dopełnieniu S, którą nazwiemy S '. W ograniczeń, można to realizować przez dodanie ograniczenie dla każdego niepusty S o rozmiarze <= k / 2 (gdzie K jest liczbą wysepek) sum_{b in S} sum_{c in S'} l_bc >= 1.

W przypadku problemu z wyspami K, kwadratami wody W i określoną maksymalną długością ścieżki N, jest to model programowania z mieszanymi liczbami całkowitymi ze O(K^2WN)zmiennymi i O(K^2WN + 2^K)ograniczeniami. Oczywiście stanie się to nie do rozwiązania, gdy rozmiar problemu stanie się duży, ale może być rozwiązany w przypadku rozmiarów, na których Ci zależy. Aby uzyskać poczucie skalowalności, zaimplementuję to w Pythonie przy użyciu pakietu pulp. Zacznijmy najpierw od mniejszej mapy 7 x 9 z 3 wyspami u dołu pytania:

import itertools
import pulp
water = {(0, 2): 2.0, (0, 3): 1.0, (0, 4): 1.0, (0, 5): 1.0, (0, 6): 2.0,
         (1, 0): 2.0, (1, 1): 9.0, (1, 2): 1.0, (1, 3): 9.0, (1, 4): 9.0,
         (1, 5): 9.0, (1, 6): 1.0, (1, 7): 9.0, (1, 8): 2.0,
         (2, 0): 1.0, (2, 1): 9.0, (2, 2): 9.0, (2, 3): 1.0, (2, 4): 9.0,
         (2, 5): 1.0, (2, 6): 9.0, (2, 7): 9.0, (2, 8): 1.0,
         (3, 0): 9.0, (3, 1): 1.0, (3, 2): 9.0, (3, 3): 9.0, (3, 4): 5.0,
         (3, 5): 9.0, (3, 6): 9.0, (3, 7): 1.0, (3, 8): 9.0,
         (4, 0): 9.0, (4, 1): 9.0, (4, 2): 1.0, (4, 3): 9.0, (4, 4): 1.0,
         (4, 5): 9.0, (4, 6): 1.0, (4, 7): 9.0, (4, 8): 9.0,
         (5, 0): 9.0, (5, 1): 9.0, (5, 2): 9.0, (5, 3): 2.0, (5, 4): 1.0,
         (5, 5): 2.0, (5, 6): 9.0, (5, 7): 9.0, (5, 8): 9.0,
         (6, 0): 9.0, (6, 1): 9.0, (6, 2): 9.0, (6, 6): 9.0, (6, 7): 9.0,
         (6, 8): 9.0}
islands = {0: [(0, 0), (0, 1)], 1: [(0, 7), (0, 8)], 2: [(6, 3), (6, 4), (6, 5)]}
N = 6

# Island borders
iborders = {}
for k in islands:
    iborders[k] = {}
    for i, j in islands[k]:
        for dx in [-1, 0, 1]:
            for dy in [-1, 0, 1]:
                if (i+dx, j+dy) in water:
                    iborders[k][(i+dx, j+dy)] = True

# Create models with specified variables
x = pulp.LpVariable.dicts("x", water.keys(), lowBound=0, upBound=1, cat=pulp.LpInteger)
pairs = [(b, c) for b in islands for c in islands if b < c]
yvals = []
for i, j in water:
    for b, c in pairs:
        for n in range(N):
            yvals.append((i, j, b, c, n))

y = pulp.LpVariable.dicts("y", yvals, lowBound=0, upBound=1)
l = pulp.LpVariable.dicts("l", pairs, lowBound=0, upBound=1)
mod = pulp.LpProblem("Islands", pulp.LpMinimize)

# Objective
mod += sum([water[k] * x[k] for k in water])

# Valid y
for k in yvals:
    i, j, b, c, n = k
    mod += y[k] <= x[(i, j)]
    if n == 0 and not (i, j) in iborders[b]:
        mod += y[k] == 0
    elif n > 0:
        mod += y[k] <= sum([y[(i+dx, j+dy, b, c, n-1)] for dx in [-1, 0, 1] for dy in [-1, 0, 1] if (i+dx, j+dy) in water])

# Valid l
for b, c in pairs:
    mod += l[(b, c)] <= sum([y[(i, j, B, C, n)] for i, j, B, C, n in yvals if (i, j) in iborders[c] and B==b and C==c])

# All islands connected (directly or indirectly)
ikeys = islands.keys()
for size in range(1, len(ikeys)/2+1):
    for S in itertools.combinations(ikeys, size):
        thisSubset = {m: True for m in S}
        Sprime = [m for m in ikeys if not m in thisSubset]
        mod += sum([l[(min(b, c), max(b, c))] for b in S for c in Sprime]) >= 1

# Solve and output
mod.solve()
for row in range(min([m[0] for m in water]), max([m[0] for m in water])+1):
    for col in range(min([m[1] for m in water]), max([m[1] for m in water])+1):
        if (row, col) in water:
            if x[(row, col)].value() > 0.999:
                print "B",
            else:
                print "-",
        else:
            print "I",
    print ""

Trwa to 1,4 sekundy do uruchomienia przy użyciu domyślnego solwera z pakietu pulp (solver CBC) i wyświetla poprawne rozwiązanie:

I I - - - - - I I 
- - B - - - B - - 
- - - B - B - - - 
- - - - B - - - - 
- - - - B - - - - 
- - - - B - - - - 
- - - I I I - - - 

Następnie rozważ pełny problem u góry pytania, czyli siatkę 13 x 14 z 7 wyspami:

water = {(i, j): 1.0 for i in range(13) for j in range(14)}
islands = {0: [(0, 0), (0, 1), (1, 0), (1, 1), (2, 0), (2, 1)],
           1: [(9, 0), (9, 1), (10, 0), (10, 1), (10, 2), (11, 0), (11, 1),
               (11, 2), (12, 0)],
           2: [(0, 7), (0, 8), (1, 7), (1, 8), (2, 7)],
           3: [(7, 7), (8, 6), (8, 7), (8, 8), (9, 7)],
           4: [(0, 11), (0, 12), (0, 13), (1, 12)],
           5: [(4, 10), (4, 11), (5, 10), (5, 11)],
           6: [(11, 8), (11, 9), (11, 13), (12, 8), (12, 9), (12, 10), (12, 11),
               (12, 12), (12, 13)]}
for k in islands:
    for i, j in islands[k]:
        del water[(i, j)]

for i, j in [(10, 7), (10, 8), (10, 9), (10, 10), (10, 11), (10, 12),
             (11, 7), (12, 7)]:
    water[(i, j)] = 20.0

N = 7

Osoby rozwiązujące MIP często uzyskują dobre rozwiązania stosunkowo szybko, a następnie spędzają dużo czasu, próbując udowodnić optymalność rozwiązania. Używając tego samego kodu solwera, co powyżej, program nie kończy się w ciągu 30 minut. Możesz jednak podać czas rozwiązującemu, aby uzyskać przybliżone rozwiązanie:

mod.solve(pulp.solvers.PULP_CBC_CMD(maxSeconds=120))

Daje to rozwiązanie o obiektywnej wartości 17:

I I - - - - - I I - - I I I 
I I - - - - - I I - - - I - 
I I - - - - - I - B - B - - 
- - B - - - B - - - B - - - 
- - - B - B - - - - I I - - 
- - - - B - - - - - I I - - 
- - - - - B - - - - - B - - 
- - - - - B - I - - - - B - 
- - - - B - I I I - - B - - 
I I - B - - - I - - - - B - 
I I I - - - - - - - - - - B 
I I I - - - - - I I - - - I 
I - - - - - - - I I I I I I 

Aby poprawić jakość otrzymywanych rozwiązań, możesz skorzystać z komercyjnego solwera MIP (jest on bezpłatny, jeśli jesteś w instytucji akademickiej i prawdopodobnie nie jest bezpłatny). Na przykład, oto wydajność Gurobi 6.0.4, ponownie z 2-minutowym limitem czasu (chociaż z dziennika rozwiązań czytamy, że solver znalazł aktualnie najlepsze rozwiązanie w ciągu 7 sekund):

mod.solve(pulp.solvers.GUROBI(timeLimit=120))

W rzeczywistości znajduje to rozwiązanie o wartości celu 16, lepsze niż to, które PO był w stanie znaleźć ręcznie!

I I - - - - - I I - - I I I 
I I - - - - - I I - - - I - 
I I - - - - - I - B - B - - 
- - B - - - - - - - B - - - 
- - - B - - - - - - I I - - 
- - - - B - - - - - I I - - 
- - - - - B - - B B - - - - 
- - - - - B - I - - B - - - 
- - - - B - I I I - - B - - 
I I - B - - - I - - - - B - 
I I I - - - - - - - - - - B 
I I I - - - - - I I - - - I 
I - - - - - - - I I I I I I 

Podejście siłowe, w pseudokodzie:

start with a horrible "best" answer
given an nxm map,
    try all 2^(n*m) combinations of bridge/no-bridge for each cell
        if the result is connected, and better than previous best, store it

return best

W C ++ można to zapisać jako

// map = linearized map; map[x*n + y] is the equivalent of map2d[y][x]
// nm = n*m
// bridged = true if bridge there, false if not. Also linearized
// nBridged = depth of recursion (= current bridge being considered)
// cost = total cost of bridges in 'bridged'
// best, bestCost = best answer so far. Initialized to "horrible"
void findBestBridges(char map[], int nm,
   bool bridged[], int nBridged, int cost, bool best[], int &bestCost) {
   if (nBridged == nm) {
      if (connected(map, nm, bridged) && cost < bestCost) {
          memcpy(best, bridged, nBridged);
          bestCost = best;
      }
      return;
   }
   if (map[nBridged] != 0) {
      // try with a bridge there
      bridged[nBridged] = true;
      cost += map[nBridged];

      // see how it turns out
      findBestBridges(map, nm, bridged, nBridged+1, cost, best, bestCost);         

      // remove bridge for further recursion
      bridged[nBridged] = false;
      cost -= map[nBridged];
   }
   // and try without a bridge there
   findBestBridges(map, nm, bridged, nBridged+1, cost, best, bestCost);
}

Po wykonaniu pierwszego połączenia (zakładam, że przekształcasz swoje mapy 2d w tablice 1d w celu ułatwienia kopiowania), bestCostbędzie zawierał koszt najlepszej odpowiedzi i bestbędzie zawierał wzór mostów, który ją daje. Jest to jednak niezwykle powolne.

Optymalizacje:

• Używając „limitu mostów” i uruchamiając algorytm zwiększania maksymalnej liczby mostów, możesz znaleźć minimalne odpowiedzi bez eksploracji całego drzewa. Znalezienie 1-mostkowej odpowiedzi, gdyby istniała, oznaczałoby O (nm) zamiast O (2 ^ nm) - drastyczna poprawa.
• Możesz uniknąć wyszukiwania (zatrzymując rekursję; jest to również nazywane „przycinaniem”) po przekroczeniu bestCost , ponieważ nie ma sensu szukać później. Jeśli nie może być lepiej, nie kop dalej.
• Powyższe przycinanie działa lepiej, jeśli spojrzysz na „dobrych” kandydatów, zanim spojrzysz na „złych” (w rzeczywistości wszystkie komórki są przeglądane w kolejności od lewej do prawej, od góry do dołu). Dobrą heurystyką byłoby traktowanie komórek, które są blisko kilku niepołączonych komponentów, jako komórek o wyższym priorytecie niż komórek, które nie są. Jednak po dodaniu heurystyki wyszukiwanie zacznie przypominać A * (i potrzebujesz też jakiejś kolejki priorytetów).
• Należy unikać powielania mostów i mostów prowadzących donikąd. Każdy most, który nie rozłącza sieci wyspowej po usunięciu, jest nadmiarowy.

Ogólny algorytm wyszukiwania, taki jak A *, umożliwia znacznie szybsze wyszukiwanie, chociaż znalezienie lepszej heurystyki nie jest prostym zadaniem. Aby uzyskać bardziej specyficzne podejście, należy wykorzystać istniejące wyniki na drzewach Steinera , zgodnie z sugestią @Gassa. Należy jednak pamiętać, że zgodnie z artykułem Gareya i Johnsona problem budowania drzew Steinera na siatkach ortogonalnych jest NP-Complete .

Jeśli wystarczy „wystarczająco dobry”, algorytm genetyczny prawdopodobnie może szybko znaleźć akceptowalne rozwiązania, o ile dodasz kilka kluczowych heurystyk dotyczących preferowanego umieszczania mostów.

Ten problem jest wariantem drzewa Steinera zwanego drzewem Steinera ważonym węzłami , specjalizującym się w określonej klasie grafów. W skrócie, drzewo Steinera ważone węzłem, biorąc pod uwagę niekierunkowy wykres ważony węzłem, w którym niektóre węzły są terminalami, znajduje najtańszy zestaw węzłów, w tym wszystkie terminale, które indukują podłączony podgraf. Niestety, w niektórych pobieżnych poszukiwaniach nie mogę znaleźć żadnego rozwiązania.

Aby sformułować jako program całkowity, utwórz zmienną 0-1 dla każdego węzła nieterminalnego, a następnie dla wszystkich podzbiorów węzłów nieterminalnych, których usunięcie z grafu początkowego rozłącza dwa terminale, wymagaj, aby suma zmiennych w podzbiorze była równa najmniej 1. To wywołuje zbyt wiele ograniczeń, więc będziesz musiał wymuszać je leniwie, używając wydajnego algorytmu łączenia węzłów (w zasadzie maksymalny przepływ), aby wykryć maksymalnie naruszone ograniczenie. Przepraszamy za brak szczegółów, ale będzie to trudne do wdrożenia, nawet jeśli znasz już programowanie liczb całkowitych.

Biorąc pod uwagę, że problem ten występuje w siatce i masz dobrze zdefiniowane parametry, podszedłbym do problemu z systematyczną eliminacją przestrzeni problemowej poprzez utworzenie minimalnego drzewa rozpinającego. Robiąc to, ma dla mnie sens, jeśli podejdziesz do tego problemu za pomocą algorytmu Prima.

Niestety, teraz napotkasz problem wyodrębnienia siatki w celu utworzenia zestawu węzłów i krawędzi… ergo prawdziwym problemem tego postu jest to, jak przekonwertować moją siatkę nxm na {V} i {E}?

Ten proces konwersji jest na pierwszy rzut oka prawdopodobnie NP-trudny ze względu na ogromną liczbę możliwych kombinacji (załóżmy, że wszystkie koszty dróg wodnych są identyczne). Aby poradzić sobie z przypadkami nakładania się ścieżek, należy rozważyć utworzenie wirtualnej wyspy.

Gdy to zrobisz, uruchom algorytm Prim i powinieneś znaleźć optymalne rozwiązanie.

Nie wierzę, że programowanie dynamiczne można tutaj skutecznie uruchomić, ponieważ nie ma zauważalnej zasady optymalności. Jeśli znajdziemy minimalny koszt między dwiema wyspami, nie musi to oznaczać, że możemy znaleźć minimalny koszt między tymi dwiema a trzecią wyspą lub inny podzbiór wysp, który będzie (zgodnie z moją definicją, aby znaleźć MST przez Prim) połączony.

Jeśli chcesz, aby kod (pseudo lub inny) przekształcił twoją siatkę w zestaw {V} i {E}, wyślij mi prywatną wiadomość, a ja zajmę się łączeniem razem implementacji.