跳转至

費用流

在看這篇文章前請先看 網絡流簡介 這篇 wiki 的定義部分。

費用流

給定一個網絡 \(G=(V,E)\)，每條邊除了有容量限制 \(c(u,v)\)，還有一個單位流量的費用 \(w(u,v)\)。

當 \((u,v)\) 的流量為 \(f(u,v)\) 時，需要花費 \(f(u,v)\times w(u,v)\) 的費用。

\(w\) 也滿足斜對稱性，即 \(w(u,v)=-w(v,u)\)。

則該網絡中總花費最小的最大流稱為 最小費用最大流，即在最大化 \(\sum_{(s,v)\in E}f(s,v)\) 的前提下最小化 \(\sum_{(u,v)\in E}f(u,v)\times w(u,v)\)。

SSP 算法

SSP（Successive Shortest Path）算法是一個貪心的算法。它的思路是每次尋找單位費用最小的增廣路進行增廣，直到圖上不存在增廣路為止。

如果圖上存在單位費用為負的圈，SSP 算法無法正確求出該網絡的最小費用最大流。此時需要先使用消圈算法消去圖上的負圈。

證明

我們考慮使用數學歸納法和反證法來證明 SSP 算法的正確性。

設流量為 \(i\) 的時候最小費用為 \(f_i\)。我們假設最初的網絡上 沒有負圈，這種情況下 \(f_0=0\)。

假設用 SSP 算法求出的 \(f_i\) 是最小費用，我們在 \(f_i\) 的基礎上，找到一條最短的增廣路，從而求出 \(f_{i+1}\)。這時 \(f_{i+1}-f_i\) 是這條最短增廣路的長度。

假設存在更小的 \(f_{i+1}\)，設它為 \(f'_{i+1}\)。因為 \(f_{i+1}-f_i\) 已經是最短增廣路了，所以 \(f'_{i+1}-f_i\) 一定對應一個經過 至少一個負圈 的增廣路。

這時候矛盾就出現了：既然存在一條經過至少一個負圈的增廣路，那麼 \(f_i\) 就不是最小費用了。因為只要給這個負圈添加流量，就可以在不增加 \(s\) 流出的流量的前提下，使 \(f_i\) 對應的費用更小。

綜上，SSP 算法可以正確求出無負圈網絡的最小費用最大流。

時間複雜度

如果使用 Bellman–Ford 算法 求解最短路，每次找增廣路的時間複雜度為 \(O(nm)\)。設該網絡的最大流為 \(f\)，則最壞時間複雜度為 \(O(nmf)\)。事實上，SSP 算法是 偽多項式時間 的。

為什麼 SSP 算法是偽多項式時間的？

SSP 算法的時間複雜度有 \(O(nmf)\) 的上界，這是一個關於值域的多項式，所以是偽多項式時間的。

可以構造 \(m=n^2,f=2^{n/2}\) 的網絡1使得 SSP 算法的時間複雜度達到 \(O(n^3 2^{n/2})\)，所以 SSP 算法不是多項式時間的。

實現

只需將 EK 算法或 Dinic 算法中找增廣路的過程，替換為用最短路算法尋找單位費用最小的增廣路即可。

基於 EK 算法的實現 
 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] = 0;
    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] = 1;
    }
  }
  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;
  }
}

// 調用：while(spfa())update();
基於 Dinic 算法的實現 
 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>

const int N = 5e3 + 5, M = 1e5 + 5;
const 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] = 1;
  while (!q.empty()) {
    int u = q.front();
    q.pop(), vis[u] = 0;
    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] = 1;
      }
    }
  }
  return dis[t] != INF;
}

int dfs(int u, int t, int flow) {
  if (u == t) return flow;
  vis[u] = 1;
  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] = 0;
  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 原始對偶算法

用 Bellman–Ford 求解最短路的時間複雜度為 \(O(nm)\)，無論在稀疏圖上還是稠密圖上都不及 Dijkstra 算法2。但網絡上存在單位費用為負的邊，因此無法直接使用 Dijkstra 算法。

Primal-Dual 原始對偶算法的思路與 Johnson 全源最短路徑算法 類似，通過為每個點設置一個勢能，將網絡上所有邊的費用（下面簡稱為邊權）全部變為非負值，從而可以應用 Dijkstra 算法找出網絡上單位費用最小的增廣路。

首先跑一次最短路，求出源點到每個點的最短距離（也是該點的初始勢能）\(h_i\)。接下來和 Johnson 算法一樣，對於一條從 \(u\) 到 \(v\)，單位費用為 \(w\) 的邊，將其邊權重置為 \(w+h_u-h_v\)。

可以發現，這樣設置勢能後新網絡上的最短路徑和原網絡上的最短路徑一定對應。證明在介紹 Johnson 算法時已經給出，這裏不再展開。

與常規的最短路問題不同的是，每次增廣後圖的形態會發生變化，這種情況下各點的勢能需要更新。

如何更新呢？先給出結論，設增廣後從源點到 \(i\) 號點的最短距離為 \(d'_i\)（這裏的距離為重置每條邊邊權後得到的距離），只需給 \(h_i\) 加上 \(d'_i\) 即可。下面我們證明，這樣更新邊權後，圖上所有邊的邊權均為非負。

容易發現，在一輪增廣後，由於一些 \((i,j)\) 邊在增廣路上，殘量網絡上會相應多出一些 \((j,i)\) 邊，且一定會滿足 \(d'_i+(w(i,j)+h_i-h_j)=d'_j\)（否則 \((i,j)\) 邊就不會在增廣路上了）。稍作變形後可以得到 \(w(j,i)+(h_j+d'_j)-(h_i+d'_i)=0\)。因此新增的邊的邊權非負。

而對於原有的邊，在增廣前，\(d'_i+(w(i,j)+h_i-h_j) - d'_j \geq 0\)，因此 \(w(i,j)+(d'_i+h_i)-(d'_j+h_j) \geq 0\)，即用 \(h_i+d'_i\) 作為新勢能並不會使 \((i,j)\) 的邊權變為負。

綜上，增廣後所有邊的邊權均非負，使用 Dijkstra 算法可以正確求出圖上的最短路。

參考代碼 
  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>
#define 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();  // 先求出初始勢能
  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;
}

習題

參考資料與註釋

  1. 詳細構造方法可以參考 min_25 的博客。 ↩

  2. 在稀疏圖上使用堆優化可以做到 \(O(m \log n)\) 的時間複雜度，而在稠密圖上不使用堆優化，可以做到 \(O(n^2)\) 的時間複雜度。 ↩

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