Min Cost Max Flow algorithm that focuses on equal distribution of flow across all edges, as much as possible

Question

My use case requires solving min cost max flow problem. I am looking for an algorithm that can satisfy the following restriction. I want to add a special restriction to finding the min cost solution. The restriction is that the cost should be calculated based on the square of the flow running through edge, not a unit cost. This restriction will force the algorithm to distribute the flow more equally.

Thank you.

Answer 1

I did two implementations in Python: one using OR-Tools and line search, the other using cvxpy. The former should be easy to transliterate to C++ or another OR-Tools language, but it's an approximation.

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)