Was sind die Mindestkosten für die Verbindung aller Inseln?

Es ist ein Gitter mit einer Größe N x M . Einige Zellen sind Inseln, die mit '0' bezeichnet sind, und die anderen sind Wasser . Auf jeder Wasserzelle befindet sich eine Nummer, die die Kosten einer auf dieser Zelle hergestellten Brücke angibt. Sie müssen die Mindestkosten finden, für die alle Inseln verbunden werden können. Eine Zelle ist mit einer anderen Zelle verbunden, wenn sie eine Kante oder einen Scheitelpunkt teilt.

Mit welchem ​​Algorithmus kann dieses Problem gelöst werden? Was kann als Brute-Force-Ansatz verwendet werden, wenn die Werte von N, M sehr klein sind, beispielsweise NxM <= 100?

Beispiel : In dem gegebenen Bild zeigen grüne Zellen Inseln an, blaue Zellen Wasser und hellblaue Zellen die Zellen, auf denen eine Brücke gebaut werden soll. Für das folgende Bild lautet die Antwort also 17 .

Anfangs dachte ich daran, alle Inseln als Knoten zu markieren und jedes Inselpaar durch eine kürzeste Brücke zu verbinden. Dann könnte das Problem auf Minimum Spanning Tree reduziert werden, aber bei diesem Ansatz habe ich den Fall übersehen, in dem sich Kanten überlappen. In der folgenden Abbildung beträgt der kürzeste Abstand zwischen zwei beliebigen Inseln 7 (gelb markiert). Wenn Sie also Minimum Spanning Trees verwenden, lautet die Antwort 14 , die Antwort sollte jedoch lauten 11 (hellblau markiert) sein.

Um dieses Problem anzugehen, würde ich ein ganzzahliges Programmierframework verwenden und drei Sätze von Entscheidungsvariablen definieren:

• x_ij : Eine binäre Indikatorvariable dafür, ob wir am Wasserstandort eine Brücke bauen (i, j).
• y_ijbcn : Ein binärer Indikator dafür, ob der Wasserort (i, j) der n-te Ort ist, der Insel b mit Insel c verbindet.
• l_bc : Eine binäre Indikatorvariable für die direkte Verknüpfung der Inseln b und c (auch bekannt als Sie können nur auf Brückenquadraten von b nach c gehen).

Für die Brückenbaukosten c_ij ist der zu minimierende Zielwertsum_ij c_ij * x_ij . Wir müssen dem Modell die folgenden Einschränkungen hinzufügen:

• Wir müssen sicherstellen, dass die Variablen y_ijbcn gültig sind. Wir können immer nur einen Wasserplatz erreichen, wenn wir dort eine Brücke bauen, also y_ijbcn <= x_ijfür jeden Wasserstandort (i, j). Ferner y_ijbc1muss 0 gleich sein, wenn (i, j) nicht an Insel b grenzt. Schließlich kann für n> 1 y_ijbcnnur verwendet werden, wenn in Schritt n-1 ein benachbarter Wasserort verwendet wurde. Definieren N(i, j)das Wasser Quadrate zu benachbarten (i, j), dann ist dies äquivalent zuy_ijbcn <= sum_{(l, m) in N(i, j)} y_lmbc(n-1) .
• Wir müssen sicherstellen, dass die Variablen l_bc nur gesetzt werden, wenn b und c verknüpft sind. Wenn wir I(c)die Orte definieren , die an Insel c grenzen, kann dies mit erreicht werdenl_bc <= sum_{(i, j) in I(c), n} y_ijbcn .
• Wir müssen sicherstellen, dass alle Inseln direkt oder indirekt miteinander verbunden sind. Dies kann auf folgende Weise erreicht werden: Für jede nicht leere richtige Teilmenge S von Inseln muss mindestens eine Insel in S mit mindestens einer Insel im Komplement von S verbunden sein, die wir S 'nennen. In Einschränkungen können wir dies implementieren, indem wir eine Einschränkung für jede nicht leere Menge S der Größe <= K / 2 hinzufügen (wobei K die Anzahl der Inseln ist) sum_{b in S} sum_{c in S'} l_bc >= 1.

Für eine Probleminstanz mit K Inseln, W Wasserquadraten und der angegebenen maximalen Pfadlänge N ist dies ein gemischtes ganzzahliges Programmiermodell mit O(K^2WN)Variablen und O(K^2WN + 2^K)Einschränkungen. Offensichtlich wird dies unlösbar, wenn die Problemgröße groß wird, aber es kann für die Größen, die Sie interessieren, lösbar sein. Um ein Gefühl für die Skalierbarkeit zu bekommen, werde ich es mit dem Pulp-Paket in Python implementieren. Beginnen wir zunächst mit der kleineren 7 x 9-Karte mit 3 Inseln am Ende der Frage:

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 ""

Die Ausführung mit dem Standardlöser aus dem Zellstoffpaket (dem CBC-Löser) dauert 1,4 Sekunden und gibt die richtige Lösung aus:

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

Betrachten Sie als nächstes das vollständige Problem oben in der Frage, nämlich ein 13 x 14-Raster mit 7 Inseln:

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

MIP-Löser erhalten oft relativ schnell gute Lösungen und verbringen dann viel Zeit damit, die Optimalität der Lösung zu beweisen. Unter Verwendung des gleichen Solver-Codes wie oben wird das Programm nicht innerhalb von 30 Minuten abgeschlossen. Sie können dem Solver jedoch eine Zeitüberschreitung gewähren, um eine ungefähre Lösung zu erhalten:

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

Dies ergibt eine Lösung mit dem Zielwert 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 

Um die Qualität der Lösungen zu verbessern, die Sie erhalten, können Sie einen kommerziellen MIP-Solver verwenden (dies ist kostenlos, wenn Sie an einer akademischen Einrichtung sind und ansonsten wahrscheinlich nicht kostenlos). Hier ist zum Beispiel die Leistung von Gurobi 6.0.4, ebenfalls mit einem Zeitlimit von 2 Minuten (obwohl wir aus dem Lösungsprotokoll gelesen haben, dass der Löser innerhalb von 7 Sekunden die derzeit beste Lösung gefunden hat):

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

Dies findet tatsächlich eine Lösung mit dem objektiven Wert 16, eine, die besser ist, als das OP von Hand finden konnte!

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 

Ein Brute-Force-Ansatz im Pseudocode:

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

In C ++ könnte dies wie folgt geschrieben werden

// 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);
}

Nach einem ersten Anruf (ich gehe davon aus, dass Sie Ihre 2D-Karten in 1D-Arrays umwandeln, um das Kopieren zu vereinfachen), bestCost werden die Kosten für die beste Antwort und bestdas Muster der Brücken, die sie ergeben, enthalten sein. Dies ist jedoch extrem langsam.

Optimierungen:

• Wenn Sie ein "Brückenlimit" verwenden und den Algorithmus zum Erhöhen der maximalen Anzahl von Brücken ausführen, können Sie minimale Antworten finden, ohne den gesamten Baum zu erkunden. Eine 1-Brücken-Antwort zu finden, falls vorhanden, wäre O (nm) anstelle von O (2 ^ nm) - eine drastische Verbesserung.
• Sie können die Suche vermeiden (indem Sie die Rekursion stoppen; dies wird auch als "Beschneiden" bezeichnet), sobald Sie überschritten haben bestCost, da es keinen Sinn macht, danach weiter zu suchen. Wenn es nicht besser werden kann, graben Sie nicht weiter.
• Das obige Beschneiden funktioniert besser, wenn Sie sich "gute" Kandidaten ansehen, bevor Sie sich "schlechte" ansehen (so wie es ist, werden alle Zellen von links nach rechts und von oben nach unten betrachtet). Eine gute Heuristik wäre, Zellen, die sich in der Nähe mehrerer nicht verbundener Komponenten befinden, als höher priorisiert zu betrachten als Zellen, die dies nicht sind. Sobald Sie jedoch Heuristiken hinzufügen, ähnelt Ihre Suche A * (und Sie benötigen auch eine Art Prioritätswarteschlange).
• Doppelte Brücken und Brücken ins Nirgendwo sind zu vermeiden. Jede Brücke, die das Inselnetzwerk nicht trennt, wenn sie entfernt wird, ist redundant.

Ein allgemeiner Suchalgorithmus wie A * ermöglicht eine viel schnellere Suche, obwohl das Finden besserer Heuristiken keine einfache Aufgabe ist. Für einen problemspezifischeren Ansatz ist die Verwendung vorhandener Ergebnisse für Steiner-Bäume , wie von @Gassa vorgeschlagen, der richtige Weg. Beachten Sie jedoch, dass das Problem beim Bauen von Steiner-Bäumen auf orthogonalen Gittern laut diesem Artikel von Garey und Johnson NP-Complete ist .

Wenn "gut genug" ausreicht, kann ein genetischer Algorithmus wahrscheinlich schnell akzeptable Lösungen finden, solange Sie einige wichtige Heuristiken für die bevorzugte Platzierung der Brücke hinzufügen.

Dieses Problem ist eine Variante des Steiner-Baums, der als knotengewichteter Steiner-Baum bezeichnet wird und auf eine bestimmte Klasse von Graphen spezialisiert ist. Kompakterweise findet der knotengewichtete Steiner-Baum bei einem knotengewichteten ungerichteten Graphen, bei dem einige Knoten Terminals sind, den billigsten Satz von Knoten, einschließlich aller Terminals, die einen verbundenen Untergraphen induzieren. Leider kann ich bei einigen flüchtigen Suchen keine Löser finden.

Um als ganzzahliges Programm zu formulieren, erstellen Sie für jeden nicht-terminalen Knoten eine 0-1-Variable. Für alle Teilmengen von nicht-terminalen Knoten, deren Entfernung aus dem Startdiagramm zwei Terminals trennt, muss die Summe der Variablen in der Teilmenge bei sein mindestens 1. Dies führt zu viel zu vielen Einschränkungen, sodass Sie diese träge durchsetzen müssen, indem Sie einen effizienten Algorithmus für die Knotenkonnektivität (im Grunde maximaler Fluss) verwenden, um eine maximal verletzte Einschränkung zu erkennen. Entschuldigen Sie den Mangel an Details, aber die Implementierung wird schwierig, selbst wenn Sie bereits mit der Ganzzahlprogrammierung vertraut sind.

Angesichts der Tatsache, dass dieses Problem in einem Raster auftritt und Sie genau definierte Parameter haben, würde ich mich dem Problem mit der systematischen Beseitigung des Problemraums durch Erstellen eines minimalen Spannbaums nähern. Dabei ist es für mich sinnvoll, dieses Problem mit dem Prim-Algorithmus anzugehen.

Leider haben Sie jetzt das Problem, das Raster zu abstrahieren, um eine Reihe von Knoten und Kanten zu erstellen. Das eigentliche Problem dieses Beitrags ist also, wie ich mein nxm-Raster in {V} und {E} konvertiere.

Dieser Umwandlungsprozess ist auf den ersten Blick aufgrund der Vielzahl möglicher Kombinationen wahrscheinlich NP-hart (vorausgesetzt, alle Wasserstraßenkosten sind identisch). Um Fälle zu behandeln, in denen sich Pfade überschneiden, sollten Sie eine virtuelle Insel erstellen.

Wenn dies erledigt ist, führen Sie den Prim-Algorithmus aus und Sie sollten zur optimalen Lösung gelangen.

Ich glaube nicht, dass die dynamische Programmierung hier effektiv ausgeführt werden kann, da es kein beobachtbares Prinzip der Optimalität gibt. Wenn wir die Mindestkosten zwischen zwei Inseln finden, bedeutet dies nicht unbedingt, dass wir die Mindestkosten zwischen diesen beiden und der dritten Insel oder einer anderen Teilmenge von Inseln finden können (nach meiner Definition, um die MST über Prim zu finden). in Verbindung gebracht.

Wenn Sie möchten, dass Code (Pseudo oder auf andere Weise) Ihr Raster in eine Menge von {V} und {E} konvertiert, senden Sie mir bitte eine private Nachricht, und ich werde eine Implementierung zusammenfügen.