费用流 Dijkstra 原始对偶方法(primal-dual method)

简单叙述用Dijkstra求费用流

Dijkstra不能求有负权边的最短路。
类似于Johnson算法，我们也可以设计一个势函数，以满足在与原图等价的新图中的边权非负。

但是这个算法并不能处理有负圈的情况（可能需要消圈算法）。

对网络\(G\)中的每一个点设置一个势函数\(h(u)\)，在任意残留网络G'的任意边\((u, v)\)都需要满足\(w_{u, v} + h(u) - h(v) \ge 0\)。
令图G的对偶图（不知道能不能这么说）为\(G'\)，其对应的边\((u, v)\)的权值为\(w'_{u, v} = w_{u, v} + h(u) + h(v)\)
对于原图中的任意一条路径\((u_1, u_2, \cdots, u_k)\)，它在\(G\)中的权值为\(w_{u_1, u_2} + w_{u_2, u_3} + \cdots + w_{u_{k - 1}, u_k}\)，在\(G'\)中的权值为
\(w'_{u_1, u_2} + w'_{u_2, u_3} + \cdots + w'_{u_{k - 1}, u_k}\)
\(= w_{u_1, u_2} + w_{u_2, u_3} + \cdots + w_{u_{k - 1}, u_k} + h(u_1) - h(u_2) + h(u_2) - h(u_3) + \cdots + h(u_{k - 1}) - h(u_k)\)
\(= w_{u_1, u_2} + w_{u_2, u_3} + \cdots + w_{u_{k - 1}, u_k} + h(u_1) - h(u_k)\)

所以，我们在\(G'\)求出的路径都可以对应到\(G\)上，令\(dist_{u, v}\)为图\(G\)点\(u\)到点\(v\)的最短路径，\(dist'_{u, v}\)为图\(G'\)点\(u\)到点\(v\)的最短路径，显然有\(dist_{u, v} = dist'_{u, v} - h(u) + h(v)\)
所以我们只需要求\(G'\)的最短路径，就能对应回原图的最短路径。

若网络\(G\)初始边权非负，我们可令\(h(u) = 0\)
否则我们令\(h(u) = dist_{s, u}\)，这个可以用Bellman-Ford算法解决。（这与Johnson算法是一模一样的）

我们考虑怎么维护势函数\(h(u)\)。
令网路\(G\)上\(dist_{s, u} = d_u\)，网路\(G'\)上\(dist'_{s, u} = d'_u\)
令残余网路\(G'\)上新的势函数为\(h'(u)\)。
对于残余网络上的一条边\((u, v)\)，有两种可能：

\((u, v) \in G\)，那么有\(d_u + w_{u, v} + h(u) - h(v) \ge d_v\)
移项得\(w_{u, v} + (h(u) + d_u) - (h(v) + d_v) \ge 0\)

\((u, v) \in G\)，那么\((v, u) \in G的增广路\)，就有\(d_v + w_{v, u} + h(v) - h(u) = d_u\)
移项得\(-w_{v, u} + (h(u) + d_u) - (h(v) - d_v) = 0\)
由\(w_{u, v} = -w_{v, u}\)可得\(w_{u, v} + (h(u) + d_u) - (h(v) - d_v) = 0\)
也就有\(w_{u, v} + (h(u) + d_u) - (h(v) - d_v) \ge 0\)

所以我们不妨令\(h'(u) = h(u) + d_u\)，这就维护好了势函数\(h(u)\)

luogu P3381 【模板】最小费用最大流

单路增广

#include <bits/stdc++.h>

using namespace std;

typedef pair<int, int> PII;

#define fi first
#define se second
#define mp make_pair

const int N = 5005, M = 50005;
const int INF = ~0u >> 1;

int n, m, s, t;
int tot = -1, head[N];
struct Edge {
    int p, nxt, c, w;
    Edge(int p = 0, int nxt = 0, int c = 0, int w = 0) : p(p), nxt(nxt), c(c), w(w) {}
} edge[M * 2];
inline void Add_Edge(int u, int v, int c, int w) {
    edge[++tot] = Edge(v, head[u], c, w);
    head[u] = tot;
    return;
}

int h[N], d[N], dc[N], pr[N];
priority_queue<PII> pq;

bool Dijkstra(int s, int t, int &mf, int &mc) {
    fill(d + 1, d + n + 1, INF);
    dc[s] = INF;
    d[s] = 0;
    pr[s] = -1;
    pq.push(mp(0, s));
    while (!pq.empty()) {
        PII cur = pq.top();
        pq.pop();
        int u = cur.se;
        if (-cur.fi > d[u]) continue;
        for (int i = head[u]; ~i; i = edge[i].nxt) {
            int v = edge[i].p, c = edge[i].c, w = edge[i].w + h[u] - h[v];
            if (!c) continue;
            if (d[v] > d[u] + w) {
                d[v] = d[u] + w;
                pr[v] = i;
                dc[v] = min(dc[u], c);
                pq.push(mp(-d[v], v));
            }
        }
    }
    if (d[t] == INF) return 0;
    for (int i = 1; i <= n; ++i)
        if (d[i] < INF) h[i] += d[i];
    int c = dc[t];
    mf += c;
    mc += c * h[t];
    for (int x = t; ~pr[x]; x = edge[pr[x] ^ 1].p) {
        edge[pr[x]].c -= c;
        edge[pr[x] ^ 1].c += c;
    }
    return 1;
}

int main() {
    scanf("%d%d%d%d", &n, &m, &s, &t);
    memset(head, -1, sizeof(head));
    for (int i = 1; i <= m; ++i) {
        int u, v, c, w;
        scanf("%d%d%d%d", &u, &v, &c, &w);
        Add_Edge(u, v, c, w);
        Add_Edge(v, u, 0, -w);
    }
    int mf = 0, mc = 0;
    while (Dijkstra(s, t, mf, mc));
    printf("%d %d\n", mf, mc);
    return 0;
}

多路增广（复杂度不变）

#include <bits/stdc++.h>

using namespace std;

typedef pair<int, int> PII;

#define fi first
#define se second
#define mp make_pair

const int N = 5005, M = 50005;
const int INF = ~0u >> 1;

int n, m, s, t;
int tot = -1, head[N], cur[N];
struct Edge {
    int p, nxt, c, w;
    Edge(int p = 0, int nxt = 0, int c = 0, int w = 0) : p(p), nxt(nxt), c(c), w(w) {}
} edge[M * 2];
inline void Add_Edge(int u, int v, int c, int w) {
    edge[++tot] = Edge(v, head[u], c, w);
    head[u] = tot;
    return;
}

int h[N], d[N];
priority_queue<PII> pq;

bool Dijkstra(int s, int t) {
    fill(d + 1, d + n + 1, INF);
    pq.push(mp(d[s] = 0, s));
    while (!pq.empty()) {
        PII cur = pq.top();
        pq.pop();
        int u = cur.se;
        if (-cur.fi > d[u]) continue;
        for (int i = head[u]; ~i; i = edge[i].nxt) {
            int v = edge[i].p, c = edge[i].c, w = edge[i].w + h[u] - h[v];
            if (c && d[v] > d[u] + w) pq.push(mp(-(d[v] = d[u] + w), v));
        }
    }
    return d[t] < INF;
}

int v[N];

int DFS(int u, int c, int t) {
    if (u == t) return c;
    int r = c;
    v[u] = 1;
    for (int &i = cur[u]; ~i && r; i = edge[i].nxt) {
        int v = edge[i].p, c = edge[i].c, w = edge[i].w + h[u] - h[v];
        if (!::v[v] && c && d[u] + w == d[v]) {
            int x = DFS(v, min(r, c), t);
            r -= x;
            edge[i].c -= x;
            edge[i ^ 1].c += x;
        }
    }
    v[u] = 0;
    return c - r;
}

int main() {
    scanf("%d%d%d%d", &n, &m, &s, &t);
    memset(head, -1, sizeof(head));
    for (int i = 1; i <= m; ++i) {
        int u, v, c, w;
        scanf("%d%d%d%d", &u, &v, &c, &w);
        Add_Edge(u, v, c, w);
        Add_Edge(v, u, 0, -w);
    }
    int mf = 0, mc = 0;
    while (Dijkstra(s, t)) {
        memcpy(cur, head, sizeof(cur));
        int c = DFS(s, INF, t);
        for (int i = 1; i <= n; ++i)
            if (d[i] < INF) h[i] += d[i];
        mf += c;
        mc += c * h[t];
    }
    printf("%d %d\n", mf, mc);
    return 0;
}