连接所有岛屿的最低成本是多少？
What is the minimum cost to connect all the islands?
有一个大小为 N x M 的网格。一些单元格是 岛屿 ，用“0”表示，其他单元格是 water。每个水格上都有一个数字，表示在该格上建造桥梁的成本。您必须找到可以连接所有岛屿的最低成本。如果一个单元格共享一条边或一个顶点，则该单元格连接到另一个单元格。
用什么算法可以解决这个问题？如果 N、M 的值非常小，例如 NxM <= 100，可以使用什么作为蛮力方法？
示例：在给定的图像中，绿色单元格表示岛屿，蓝色单元格表示水，浅蓝色单元格表示应在其上建造桥梁的单元格。因此对于下图，答案将是 17.
最初我想到的是将所有岛屿标记为节点，并通过最短的桥梁将每对岛屿连接起来。然后问题可以减少到最小生成树，但在这种方法中我错过了边缘重叠的情况。 例如，在下图中，任意两个岛之间的最短距离为7（黄色标记），因此通过使用最小生成树答案是 14，但答案应该是 11（以浅蓝色标记）。
伪代码中的蛮力方法：
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
在 C++ 中，这可以写成
// 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);
}
第一次调用后（我假设您正在将 2d 地图转换为 1d 数组以便于复制），bestCost 将包含最佳答案的成本，best将包含产生它的桥梁模式。但是，这非常慢。
优化：
- 通过使用 "bridge limit" 和 运行 增加最大桥数的算法，您无需探索整棵树即可找到最少的答案。找到一个单桥答案（如果存在）将是 O(nm) 而不是 O(2^nm) - 一个巨大的改进。
- 一旦超过
bestCost，您就可以避免搜索（通过停止递归；这也称为"pruning"），因为之后继续寻找是没有意义的。如果无法改善，请不要继续挖掘。
- 如果在查看 "bad" 个候选对象之前先查看 "good" 个候选对象，上述修剪效果会更好（实际上，单元格都是从左到右、从上到-底顺序）。一个好的启发式方法是将靠近几个未连接组件的单元格视为比不靠近的单元格具有更高的优先级。但是，一旦您添加了启发式算法，您的搜索就会开始类似于 A*（并且您也需要某种优先级队列）。
- 要避免重复的桥梁和无处可去的桥梁。任何在移除后不断开孤岛网络的网桥都是多余的。
通用搜索算法（例如 A*）可以实现更快的搜索，尽管找到更好的启发式算法并不是一项简单的任务。对于更具体的问题的方法，按照@Gassa 的建议，使用 Steiner 树 上的现有结果是可行的方法。但是请注意，根据 Garey 和 Johnson 的 论文。
，在正交网格上构建斯坦纳树的问题是 NP 完全问题
如果 "good enough" 就足够了，遗传算法可能可以快速找到可接受的解决方案，只要您添加一些关于首选桥放置的关键启发式算法。
这个问题是斯坦纳树的一个变体，叫做节点加权斯坦纳树，专门针对某个class图。紧凑地，节点加权斯坦纳树是，给定一个节点加权无向图，其中一些节点是终端，找到最便宜的一组节点，包括所有终端，从而导出一个连接的子图。可悲的是，我似乎无法在一些粗略的搜索中找到任何求解器。
为了制定一个整数规划，为每个非终端节点制作一个0-1变量，然后对于从起始图中删除断开两个终端的非终端节点的所有子集，要求变量之和在子集中至少为 1。这会引入太多约束，因此您必须懒惰地强制执行它们，使用有效的节点连接算法（基本上是最大流）来检测最大违反约束。很抱歉缺少详细信息，但即使您已经熟悉整数规划，实施起来也会很痛苦。
鉴于此问题发生在网格中，并且您已明确定义参数，我将通过创建最小生成树来系统地消除问题 space 来解决问题。在这样做的过程中，如果您使用 Prim 算法解决这个问题，对我来说很有意义。
不幸的是，您现在 运行 陷入了抽象网格以创建一组节点和边的问题...因此 post 的真正问题是 如何将我的 n x m 网格转换为 {V} 和 {E}？
这个转换过程乍一看很可能是 NP-Hard，因为可能的组合数量非常多（假设所有水路成本都相同）。要处理路径重叠的情况，您应该考虑制作一个 虚拟岛。
完成后，运行 Prim 算法，您应该会得出最优解。
我不相信动态规划在这里可以 运行 有效，因为没有可观察到的最优性原则。如果我们找到两个岛之间的最小成本，那并不一定意味着我们可以找到这两个岛和第三个岛之间的最小成本，或者另一个岛屿子集（根据我的定义，通过 Prim 找到 MST）已连接。
如果您想要代码（伪代码或其他代码）将您的网格转换为一组 {V} 和 {E}，请给我发私信，我会考虑拼接实现。
为了解决这个问题，我将使用整数规划框架并定义三组决策变量：
- x_ij：一个二元指标变量，表示我们是否在水位（i，j）建桥。
- y_ijbcn: 水域位置 (i, j) 是否为连接 b 岛和 c 岛的第 n^th 个位置的二元指标。
- l_bc：一个二元指示变量，表示岛 b 和岛 c 是否直接相连（也就是从 b 到 c 只能走桥广场）。
对于桥梁建造成本 c_ij，要最小化的 objective 值为 sum_ij c_ij * x_ij。我们需要在模型中添加以下约束：
- 我们需要确保 y_ijbcn 变量有效。如果我们在那里建造一座桥，我们总是只能到达一个水广场，所以
y_ijbcn <= x_ij 对于每个水位置 (i, j)。此外，如果 (i, j) 不与岛 b 接壤，则 y_ijbc1 必须等于 0。最后，对于 n > 1，只有在步骤 n-1 中使用了邻近的水域位置时，才能使用 y_ijbcn。将N(i, j)定义为与(i,j)相邻的水域方块，等价于y_ijbcn <= sum_{(l, m) in N(i, j)} y_lmbc(n-1).
- 我们需要确保 l_bc 变量仅在 b 和 c 链接时设置。如果我们将
I(c) 定义为与岛 c 接壤的位置，则可以通过 l_bc <= sum_{(i, j) in I(c), n} y_ijbcn. 来实现
- 我们需要确保所有岛屿都直接或间接相连。这可以通过以下方式实现：对于岛屿的每个非空真子集 S，要求 S 中的至少一个岛屿与 S 的补集中的至少一个岛屿相连，我们将其称为 S'。在约束中，我们可以通过为大小 <= K/2（其中 K 是岛数）、
sum_{b in S} sum_{c in S'} l_bc >= 1. 的每个非空集 S 添加约束来实现这一点
对于具有 K 个岛屿、W 个水域和指定的最大路径长度 N 的问题实例，这是一个具有 O(K^2WN) 个变量和 O(K^2WN + 2^K) 个约束的混合整数规划模型。显然，随着问题规模变大，这将变得棘手，但对于您关心的规模，它可能是可以解决的。为了了解可伸缩性，我将使用 pulp 包在 python 中实现它。让我们先从问题底部有 3 个岛屿的较小的 7 x 9 地图开始：
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 ""
这需要 1.4 秒才能 运行 使用 pulp 包中的默认求解器（CBC 求解器）并输出正确的解决方案：
I I - - - - - I I
- - B - - - B - -
- - - B - B - - -
- - - - B - - - -
- - - - B - - - -
- - - - B - - - -
- - - I I I - - -
接下来，考虑问题顶部的完整问题，即具有 7 个岛的 13 x 14 网格：
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 求解器通常会相对较快地获得良好的解决方案，然后花费大量时间来尝试证明解决方案的最优性。使用与上面相同的求解器代码，程序不会在 30 分钟内完成。但是，您可以为求解器提供超时以获取近似解：
mod.solve(pulp.solvers.PULP_CBC_CMD(maxSeconds=120))
这会产生一个 objective 值为 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
要提高您获得的解决方案的质量，您可以使用商业 MIP 求解器（如果您在学术机构，这是免费的，否则可能不是免费的）。例如，这是 Gurobi 6.0.4 的性能，同样有 2 分钟的时间限制（尽管从解决方案日志中我们读到求解器在 7 秒内找到了当前最佳解决方案）：
mod.solve(pulp.solvers.GUROBI(timeLimit=120))
这实际上找到了 objective 值 16 的解决方案，比 OP 手动找到的解决方案更好！
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
有一个大小为 N x M 的网格。一些单元格是 岛屿 ，用“0”表示，其他单元格是 water。每个水格上都有一个数字，表示在该格上建造桥梁的成本。您必须找到可以连接所有岛屿的最低成本。如果一个单元格共享一条边或一个顶点，则该单元格连接到另一个单元格。
用什么算法可以解决这个问题？如果 N、M 的值非常小，例如 NxM <= 100，可以使用什么作为蛮力方法？
示例：在给定的图像中，绿色单元格表示岛屿，蓝色单元格表示水，浅蓝色单元格表示应在其上建造桥梁的单元格。因此对于下图，答案将是 17.
最初我想到的是将所有岛屿标记为节点，并通过最短的桥梁将每对岛屿连接起来。然后问题可以减少到最小生成树，但在这种方法中我错过了边缘重叠的情况。 例如，在下图中，任意两个岛之间的最短距离为7（黄色标记），因此通过使用最小生成树答案是 14，但答案应该是 11（以浅蓝色标记）。
伪代码中的蛮力方法：
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
在 C++ 中，这可以写成
// 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);
}
第一次调用后（我假设您正在将 2d 地图转换为 1d 数组以便于复制），bestCost 将包含最佳答案的成本，best将包含产生它的桥梁模式。但是，这非常慢。
优化：
- 通过使用 "bridge limit" 和 运行 增加最大桥数的算法，您无需探索整棵树即可找到最少的答案。找到一个单桥答案（如果存在）将是 O(nm) 而不是 O(2^nm) - 一个巨大的改进。
- 一旦超过
bestCost，您就可以避免搜索（通过停止递归；这也称为"pruning"），因为之后继续寻找是没有意义的。如果无法改善，请不要继续挖掘。 - 如果在查看 "bad" 个候选对象之前先查看 "good" 个候选对象，上述修剪效果会更好（实际上，单元格都是从左到右、从上到-底顺序）。一个好的启发式方法是将靠近几个未连接组件的单元格视为比不靠近的单元格具有更高的优先级。但是，一旦您添加了启发式算法，您的搜索就会开始类似于 A*（并且您也需要某种优先级队列）。
- 要避免重复的桥梁和无处可去的桥梁。任何在移除后不断开孤岛网络的网桥都是多余的。
通用搜索算法（例如 A*）可以实现更快的搜索，尽管找到更好的启发式算法并不是一项简单的任务。对于更具体的问题的方法，按照@Gassa 的建议，使用 Steiner 树 上的现有结果是可行的方法。但是请注意，根据 Garey 和 Johnson 的 论文。
，在正交网格上构建斯坦纳树的问题是 NP 完全问题如果 "good enough" 就足够了，遗传算法可能可以快速找到可接受的解决方案，只要您添加一些关于首选桥放置的关键启发式算法。
这个问题是斯坦纳树的一个变体，叫做节点加权斯坦纳树，专门针对某个class图。紧凑地，节点加权斯坦纳树是，给定一个节点加权无向图，其中一些节点是终端，找到最便宜的一组节点，包括所有终端，从而导出一个连接的子图。可悲的是，我似乎无法在一些粗略的搜索中找到任何求解器。
为了制定一个整数规划，为每个非终端节点制作一个0-1变量，然后对于从起始图中删除断开两个终端的非终端节点的所有子集，要求变量之和在子集中至少为 1。这会引入太多约束，因此您必须懒惰地强制执行它们，使用有效的节点连接算法（基本上是最大流）来检测最大违反约束。很抱歉缺少详细信息，但即使您已经熟悉整数规划，实施起来也会很痛苦。
鉴于此问题发生在网格中，并且您已明确定义参数，我将通过创建最小生成树来系统地消除问题 space 来解决问题。在这样做的过程中，如果您使用 Prim 算法解决这个问题，对我来说很有意义。
不幸的是，您现在 运行 陷入了抽象网格以创建一组节点和边的问题...因此 post 的真正问题是 如何将我的 n x m 网格转换为 {V} 和 {E}？
这个转换过程乍一看很可能是 NP-Hard，因为可能的组合数量非常多（假设所有水路成本都相同）。要处理路径重叠的情况，您应该考虑制作一个 虚拟岛。
完成后，运行 Prim 算法，您应该会得出最优解。
我不相信动态规划在这里可以 运行 有效，因为没有可观察到的最优性原则。如果我们找到两个岛之间的最小成本，那并不一定意味着我们可以找到这两个岛和第三个岛之间的最小成本，或者另一个岛屿子集（根据我的定义，通过 Prim 找到 MST）已连接。
如果您想要代码（伪代码或其他代码）将您的网格转换为一组 {V} 和 {E}，请给我发私信，我会考虑拼接实现。
为了解决这个问题，我将使用整数规划框架并定义三组决策变量：
- x_ij：一个二元指标变量，表示我们是否在水位（i，j）建桥。
- y_ijbcn: 水域位置 (i, j) 是否为连接 b 岛和 c 岛的第 n^th 个位置的二元指标。
- l_bc：一个二元指示变量，表示岛 b 和岛 c 是否直接相连（也就是从 b 到 c 只能走桥广场）。
对于桥梁建造成本 c_ij，要最小化的 objective 值为 sum_ij c_ij * x_ij。我们需要在模型中添加以下约束：
- 我们需要确保 y_ijbcn 变量有效。如果我们在那里建造一座桥，我们总是只能到达一个水广场，所以
y_ijbcn <= x_ij对于每个水位置 (i, j)。此外，如果 (i, j) 不与岛 b 接壤，则y_ijbc1必须等于 0。最后，对于 n > 1，只有在步骤 n-1 中使用了邻近的水域位置时，才能使用y_ijbcn。将N(i, j)定义为与(i,j)相邻的水域方块，等价于y_ijbcn <= sum_{(l, m) in N(i, j)} y_lmbc(n-1). - 我们需要确保 l_bc 变量仅在 b 和 c 链接时设置。如果我们将
I(c)定义为与岛 c 接壤的位置，则可以通过l_bc <= sum_{(i, j) in I(c), n} y_ijbcn. 来实现
- 我们需要确保所有岛屿都直接或间接相连。这可以通过以下方式实现：对于岛屿的每个非空真子集 S，要求 S 中的至少一个岛屿与 S 的补集中的至少一个岛屿相连，我们将其称为 S'。在约束中，我们可以通过为大小 <= K/2（其中 K 是岛数）、
sum_{b in S} sum_{c in S'} l_bc >= 1. 的每个非空集 S 添加约束来实现这一点
对于具有 K 个岛屿、W 个水域和指定的最大路径长度 N 的问题实例，这是一个具有 O(K^2WN) 个变量和 O(K^2WN + 2^K) 个约束的混合整数规划模型。显然，随着问题规模变大，这将变得棘手，但对于您关心的规模，它可能是可以解决的。为了了解可伸缩性，我将使用 pulp 包在 python 中实现它。让我们先从问题底部有 3 个岛屿的较小的 7 x 9 地图开始：
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 ""
这需要 1.4 秒才能 运行 使用 pulp 包中的默认求解器（CBC 求解器）并输出正确的解决方案：
I I - - - - - I I
- - B - - - B - -
- - - B - B - - -
- - - - B - - - -
- - - - B - - - -
- - - - B - - - -
- - - I I I - - -
接下来，考虑问题顶部的完整问题，即具有 7 个岛的 13 x 14 网格：
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 求解器通常会相对较快地获得良好的解决方案，然后花费大量时间来尝试证明解决方案的最优性。使用与上面相同的求解器代码，程序不会在 30 分钟内完成。但是，您可以为求解器提供超时以获取近似解：
mod.solve(pulp.solvers.PULP_CBC_CMD(maxSeconds=120))
这会产生一个 objective 值为 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
要提高您获得的解决方案的质量，您可以使用商业 MIP 求解器（如果您在学术机构，这是免费的，否则可能不是免费的）。例如，这是 Gurobi 6.0.4 的性能，同样有 2 分钟的时间限制（尽管从解决方案日志中我们读到求解器在 7 秒内找到了当前最佳解决方案）：
mod.solve(pulp.solvers.GUROBI(timeLimit=120))
这实际上找到了 objective 值 16 的解决方案，比 OP 手动找到的解决方案更好！
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