Min Cost Max Flow 算法，尽可能地关注所有边缘上的流量均等分布

Posted 2023-02-19

【中文标题】Min Cost Max Flow 算法，尽可能地关注所有边缘上的流量均等分布

【英文标题】：Min Cost Max Flow algorithm that focuses on equal distribution of flow across all edges, as much as possible

【发布时间】：2021-11-16 03:38:27

【问题描述】：

我的用例需要解决最小成本最大流量问题。我正在寻找一种可以满足以下限制的算法。 我想添加一个特殊限制来寻找最低成本解决方案。 限制是成本应该基于流过边缘的流量的平方来计算，而不是单位成本。这种限制将迫使算法更平均地分配流量。

谢谢。

【问题讨论】：

这里有易于实现与可扩展性的权衡。你的图表有多大？

假设没那么大。可能有 1000 个节点用于任务，20 个节点用于工作人员。和最大（20000）条边。我不担心可扩展性，我只想要一些可行的东西。

我明白了，谢谢@DavidEisenstat 对 cvxpy 的建议。不幸的是，我正在使用谷歌开发人员或工具来解决最大流量问题，所以我没有一种方便的方法来迭代/利用 cvxpy。

非常感谢，我正在玩他们的代码，但他们的 DCP 规则很烦人。

【参考方案1】：

我在 Python 中做了两个实现：一个使用 OR-Tools 和行搜索，另一个使用 cvxpy。前者应该很容易转写为 C++ 或其他 OR-Tools 语言，但它是一个近似值。

import ortools.graph.pywrapgraph
import random

# Set up a random network.
n = 100
average_out_degree = 20
max_arc_capacity = 10
arcs = []
for j in range(n * average_out_degree):
    while True:
        tail = random.randrange(n)
        head = random.randrange(n)
        if tail != head:
            break
    capacity = random.randrange((n // abs(tail - head)) ** 2 + 1)
    # We need a lot of excess precision in the capacity because we round.
    arcs.append((tail, head, 10 ** 4 * capacity))
source = 0
sink = n - 1

# Initialize the line search with a max flow.
max_flow = ortools.graph.pywrapgraph.SimpleMaxFlow()
arc_indices = []
for tail, head, capacity in arcs:
    arc_indices.append(max_flow.AddArcWithCapacity(tail, head, capacity))
max_flow.Solve(source, sink)
flows = []
for arc_index in arc_indices:
    flows.append(max_flow.Flow(arc_index))
amount = max_flow.OptimalFlow()

# Improve the cost of this flow using line search.
for i in range(1000):
    print("Cost is", sum(flow ** 2 for flow in flows))
    # The gradient of the L2 cost function is linear. Therefore we can compute
    # an optimal improving direction as a min-cost flow.
    min_cost_flow = ortools.graph.pywrapgraph.SimpleMinCostFlow()
    arc_indices = []
    for (tail, head, capacity), flow in zip(arcs, flows):
        # OR-Tools requires the unit cost to be an integer.
        unit_cost = round(flow)
        arc_indices.append(
            min_cost_flow.AddArcWithCapacityAndUnitCost(tail, head, capacity, unit_cost)
        )
    min_cost_flow.SetNodeSupply(source, amount)
    min_cost_flow.SetNodeSupply(sink, -amount)
    min_cost_flow.Solve()
    flows_prime = []
    for arc_index in arc_indices:
        flows_prime.append(min_cost_flow.Flow(arc_index))
    # Now we take the best solution that is a convex combination (1 - alpha) *
    # flows + alpha * flows_prime. First we express the cost as a quadratic
    # polynomial in alpha.
    a = 0
    b = 0
    c = 0
    for flow, flow_prime in zip(flows, flows_prime):
        # Expand ((1 - alpha) * flow + alpha * flow_prime) ** 2 and get the
        # coefficients below.
        a += (flow_prime - flow) ** 2
        b += 2 * (flow_prime - flow) * flow
        c += flow ** 2
    if a == 0:
        # flows and flows_prime are the same. No point in continuing.
        break
    # Since a > 0, this is where the polynomial takes its minimum.
    alpha = -b / (2 * a)
    # Clamp alpha.
    alpha = max(0, min(1, alpha))
    # Update flows.
    for j, (flow, flow_prime) in enumerate(zip(flows, flows_prime)):
        flows[j] = (1 - alpha) * flow + alpha * flow_prime


# Comparison using cvxpy.

import cvxpy

flows = cvxpy.Variable(len(arcs))
constraints = []
excesses = [0] * n
excesses[source] = amount
excesses[sink] = -amount
for (tail, head, capacity), flow in zip(arcs, flows):
    excesses[tail] -= flow
    excesses[head] += flow
    constraints.append(flow >= 0)
    constraints.append(flow <= capacity)
for excess in excesses:
    constraints.append(excess == 0)
problem = cvxpy.Problem(cvxpy.Minimize(cvxpy.sum_squares(flows)), constraints)
problem.solve()
print(problem.value)