Skip to content

Minimum Cost Flow

Before reading this article, please first read the definition section of the Network Flow Introduction wiki.

Minimum Cost Flow

Given a network \(G=(V,E)\), each edge has not only a capacity constraint \(c(u,v)\) but also a cost per unit flow \(w(u,v)\).

When the flow on \((u,v)\) is \(f(u,v)\), the cost incurred is \(f(u,v)\times w(u,v)\).

\(w\) also satisfies antisymmetry, i.e., \(w(u,v)=-w(v,u)\).

Then, the maximum flow with the minimum total cost in this network is called the minimum cost maximum flow, which minimizes \(\sum_{(u,v)\in E}f(u,v)\times w(u,v)\) under the premise of maximizing \(\sum_{(s,v)\in E}f(s,v)\).

SSP Algorithm

The SSP (Successive Shortest Path) algorithm is a greedy algorithm. Its idea is to repeatedly find the augmenting path with the minimum unit cost and augment along it until no augmenting path exists in the graph.

If there exists a cycle with negative unit cost in the graph, the SSP algorithm cannot correctly compute the minimum cost maximum flow of the network. In this case, it is necessary to first use a cycle cancellation algorithm to eliminate negative cycles in the graph.

Proof

We use mathematical induction and proof by contradiction to demonstrate the correctness of the SSP algorithm.

Let the minimum cost for flow \(i\) be \(f_i\). We assume that initially the network has no negative cycles, in which case \(f_0=0\).

Assume that \(f_i\) computed by the SSP algorithm is the minimum cost. Based on \(f_i\), we find the shortest augmenting path to compute \(f_{i+1}\). Then \(f_{i+1}-f_i\) is the length of this shortest augmenting path.

Suppose there exists a smaller \(f_{i+1}\), denoted as \(f'_{i+1}\). Since \(f_{i+1}-f_i\) is already the shortest augmenting path, \(f'_{i+1}-f_i\) must correspond to an augmenting path that passes through at least one negative cycle.

A contradiction arises here: since there exists an augmenting path passing through at least one negative cycle, \(f_i\) is not the minimum cost. Because by adding flow to this negative cycle, we can reduce the cost corresponding to \(f_i\) without increasing the flow out of \(s\).

In conclusion, the SSP algorithm can correctly compute the minimum cost maximum flow for networks without negative cycles.

Time Complexity

If the Bellman–Ford algorithm is used to find the shortest path, each augmentation takes \(O(nm)\) time. Let the maximum flow of the network be \(f\), then the worst-case time complexity is \(O(nmf)\). In fact, the SSP algorithm runs in pseudo-polynomial time.

Why is the SSP algorithm pseudo-polynomial?

The time complexity of the SSP algorithm has an upper bound of \(O(nmf)\), which is a polynomial in the value domain, hence it is pseudo-polynomial.

One can construct a network2 with \(m=n^2,f=2^{n/2}\) such that the SSP algorithm's time complexity reaches \(O(n^3 2^{n/2})\), so the SSP algorithm is not polynomial-time.

Implementation

Simply replace the process of finding augmenting paths in the EK algorithm or Dinic's algorithm with a shortest path algorithm that finds the augmenting path with the minimum unit cost.

Implementation based on the EK algorithm 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
struct qxx {
  int nex, t, v, c;
};

qxx e[M];
int h[N], cnt = 1;

void add_path(int f, int t, int v, int c) {
  e[++cnt] = qxx{h[f], t, v, c}, h[f] = cnt;
}

void add_flow(int f, int t, int v, int c) {
  add_path(f, t, v, c);
  add_path(t, f, 0, -c);
}

int dis[N], pre[N], incf[N];
bool vis[N];

bool spfa() {
  memset(dis, 0x3f, sizeof(dis));
  queue<int> q;
  q.push(s), dis[s] = 0, incf[s] = INF, incf[t] = 0;
  while (q.size()) {
    int u = q.front();
    q.pop();
    vis[u] = false;
    for (int i = h[u]; i; i = e[i].nex) {
      const int &v = e[i].t, &w = e[i].v, &c = e[i].c;
      if (!w || dis[v] <= dis[u] + c) continue;
      dis[v] = dis[u] + c, incf[v] = min(w, incf[u]), pre[v] = i;
      if (!vis[v]) q.push(v), vis[v] = true;
    }
  }
  return incf[t];
}

int maxflow, mincost;

void update() {
  maxflow += incf[t];
  for (int u = t; u != s; u = e[pre[u] ^ 1].t) {
    e[pre[u]].v -= incf[t], e[pre[u] ^ 1].v += incf[t];
    mincost += incf[t] * e[pre[u]].c;
  }
}

// Call: while(spfa())update();
Implementation based on Dinic's algorithm 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <queue>

constexpr int N = 5e3 + 5, M = 1e5 + 5;
constexpr int INF = 0x3f3f3f3f;
int n, m, tot = 1, lnk[N], cur[N], ter[M], nxt[M], cap[M], cost[M], dis[N], ret;
bool vis[N];

void add(int u, int v, int w, int c) {
  ter[++tot] = v, nxt[tot] = lnk[u], lnk[u] = tot, cap[tot] = w, cost[tot] = c;
}

void addedge(int u, int v, int w, int c) { add(u, v, w, c), add(v, u, 0, -c); }

bool spfa(int s, int t) {
  memset(dis, 0x3f, sizeof(dis));
  memcpy(cur, lnk, sizeof(lnk));
  std::queue<int> q;
  q.push(s), dis[s] = 0, vis[s] = true;
  while (!q.empty()) {
    int u = q.front();
    q.pop(), vis[u] = false;
    for (int i = lnk[u]; i; i = nxt[i]) {
      int v = ter[i];
      if (cap[i] && dis[v] > dis[u] + cost[i]) {
        dis[v] = dis[u] + cost[i];
        if (!vis[v]) q.push(v), vis[v] = true;
      }
    }
  }
  return dis[t] != INF;
}

int dfs(int u, int t, int flow) {
  if (u == t) return flow;
  vis[u] = true;
  int ans = 0;
  for (int &i = cur[u]; i && ans < flow; i = nxt[i]) {
    int v = ter[i];
    if (!vis[v] && cap[i] && dis[v] == dis[u] + cost[i]) {
      int x = dfs(v, t, std::min(cap[i], flow - ans));
      if (x) ret += x * cost[i], cap[i] -= x, cap[i ^ 1] += x, ans += x;
    }
  }
  vis[u] = false;
  return ans;
}

int mcmf(int s, int t) {
  int ans = 0;
  while (spfa(s, t)) {
    int x;
    while ((x = dfs(s, t, INF))) ans += x;
  }
  return ans;
}

int main() {
  int s, t;
  scanf("%d%d%d%d", &n, &m, &s, &t);
  while (m--) {
    int u, v, w, c;
    scanf("%d%d%d%d", &u, &v, &w, &c);
    addedge(u, v, w, c);
  }
  int ans = mcmf(s, t);
  printf("%d %d\n", ans, ret);
  return 0;
}

Primal-Dual Algorithm

Using Bellman–Ford to find the shortest path has a time complexity of \(O(nm)\), which is inferior to Dijkstra's algorithm1 on both sparse and dense graphs. However, the network may contain edges with negative unit costs, so Dijkstra's algorithm cannot be directly applied.

The idea of the Primal-Dual algorithm is similar to the Johnson's algorithm for all-pairs shortest paths. By assigning a potential to each vertex, the cost (hereafter referred to as edge weight) of all edges in the network becomes non-negative, allowing Dijkstra's algorithm to find the augmenting path with the minimum unit cost.

First, run a shortest path algorithm to find the shortest distance (also the initial potential) \(h_i\) from the source to each vertex. Then, similar to Johnson's algorithm, for an edge from \(u\) to \(v\) with unit cost \(w\), reset its edge weight to \(w+h_u-h_v\).

It can be shown that the shortest paths in the new network correspond exactly to the shortest paths in the original network. The proof is given in the introduction to Johnson's algorithm and will not be repeated here.

Unlike conventional shortest path problems, the graph's structure changes after each augmentation. In this case, the potentials of the vertices need to be updated.

How to update them? First, let the shortest distance from the source to vertex \(i\) after augmentation be \(d'_i\) (here the distance is computed after resetting edge weights). Simply add \(d'_i\) to \(h_i\). Below we prove that after updating edge weights in this way, all edge weights in the graph become non-negative.

It is easy to see that after one augmentation, because some edges \((i,j)\) lie on the augmenting path, the residual network will correspondingly have additional edges \((j,i)\), and they must satisfy \(d'_i+(w(i,j)+h_i-h_j)=d'_j\) (otherwise edge \((i,j)\) would not be on the augmenting path). After slight rearrangement, we get \(w(j,i)+(h_j+d'_j)-(h_i+d'_i)=0\). Therefore, the weights of the newly added edges are non-negative.

For the original edges, before augmentation, \(d'_i+(w(i,j)+h_i-h_j) - d'_j \geq 0\), so \(w(i,j)+(d'_i+h_i)-(d'_j+h_j) \geq 0\), meaning that using \(h_i+d'_i\) as the new potential does not make the weight of \((i,j)\) negative.

In conclusion, after augmentation, all edge weights are non-negative, and Dijkstra's algorithm can correctly find the shortest paths in the graph.

Reference Code 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <queue>
constexpr int INF = 0x3f3f3f3f;
using namespace std;

struct edge {
  int v, f, c, next;
} e[100005];

struct node {
  int v, e;
} p[10005];

struct mypair {
  int dis, id;

  bool operator<(const mypair& a) const { return dis > a.dis; }

  mypair(int d, int x) { dis = d, id = x; }
};

int head[5005], dis[5005], vis[5005], h[5005];
int n, m, s, t, cnt = 1, maxf, minc;

void addedge(int u, int v, int f, int c) {
  e[++cnt].v = v;
  e[cnt].f = f;
  e[cnt].c = c;
  e[cnt].next = head[u];
  head[u] = cnt;
}

bool dijkstra() {
  priority_queue<mypair> q;
  for (int i = 1; i <= n; i++) dis[i] = INF;
  memset(vis, 0, sizeof(vis));
  dis[s] = 0;
  q.push(mypair(0, s));
  while (!q.empty()) {
    int u = q.top().id;
    q.pop();
    if (vis[u]) continue;
    vis[u] = 1;
    for (int i = head[u]; i; i = e[i].next) {
      int v = e[i].v, nc = e[i].c + h[u] - h[v];
      if (e[i].f && dis[v] > dis[u] + nc) {
        dis[v] = dis[u] + nc;
        p[v].v = u;
        p[v].e = i;
        if (!vis[v]) q.push(mypair(dis[v], v));
      }
    }
  }
  return dis[t] != INF;
}

void spfa() {
  queue<int> q;
  memset(h, 63, sizeof(h));
  h[s] = 0, vis[s] = 1;
  q.push(s);
  while (!q.empty()) {
    int u = q.front();
    q.pop();
    vis[u] = 0;
    for (int i = head[u]; i; i = e[i].next) {
      int v = e[i].v;
      if (e[i].f && h[v] > h[u] + e[i].c) {
        h[v] = h[u] + e[i].c;
        if (!vis[v]) {
          vis[v] = 1;
          q.push(v);
        }
      }
    }
  }
}

int main() {
  scanf("%d%d%d%d", &n, &m, &s, &t);
  for (int i = 1; i <= m; i++) {
    int u, v, f, c;
    scanf("%d%d%d%d", &u, &v, &f, &c);
    addedge(u, v, f, c);
    addedge(v, u, 0, -c);
  }
  spfa();  // First compute initial potentials
  while (dijkstra()) {
    int minf = INF;
    for (int i = 1; i <= n; i++) h[i] += dis[i];
    for (int i = t; i != s; i = p[i].v) minf = min(minf, e[p[i].e].f);
    for (int i = t; i != s; i = p[i].v) {
      e[p[i].e].f -= minf;
      e[p[i].e ^ 1].f += minf;
    }
    maxf += minf;
    minc += minf * h[t];
  }
  printf("%d %d\n", maxf, minc);
  return 0;
}

Exercises

References and Notes

  1. Using heap optimization on sparse graphs can achieve \(O(m \log n)\) time complexity, while on dense graphs, without heap optimization, \(O(n^2)\) time complexity can be achieved. ↩

  2. For detailed construction methods, refer to min_25's blog. ↩

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用