费用流

对于一个网络图 \(G=(V,E)\)，每条边都被额外赋予了一个权值：费用 \(c(u,v)\)。一个流的总费用是 \(\sum_{(u,v)\in E}{f(u,v)c(u,v)}\)。

最小费用最大流

求解最小费用最大流的方法仍然是增广路算法。这个算法要求网络图初始时没有负圈。

引理

如果每次增广时都选择代价最小的增广路，那么求出的一定是最小费用最大流。

证明

考虑数学归纳法。设 \(f_i\) 表示流量为 \(i\) 时的最小费用。首先任意 \(f_i\) 一定不包含负圈。然后找最小代价的增广路进行增广显然也不会产生负圈。而且只要没有负圈就一定是最小费用的流。

spfa

根据引理，我们可以使用 Dinic+spfa 实现最小费用最大流。

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

#include<iostream>
using namespace std;
const int N = 5e3 + 10;
const int M = 5e4 + 10;
const int INF = (int)0x3f3f3f3f3f3f3f3f;

struct Edge {
    int v, w, c;
    int next;
} pool[2 * M];
int ne = 1, head[N];

void addEdge(int u, int v, int w, int c) {
    pool[++ne] = {v, w, c, head[u]};
    head[u] = ne;
}

int n, m, s, t;

namespace dinic {
    int dis[N], dep[N], vis[N], now[N];
    int que[M], hd, tl;
    int ans;
    bool spfa() {
        for(int i = 1; i <= n; i++) {
            vis[i] = 0;
            dis[i] = INF;
            now[i] = head[i];
        }
        hd = 1, tl = 0;
        que[++tl] = s;
        dis[s] = 0;
        while(hd <= tl) {
            int u = que[hd++];
            vis[u] = 0;
            for(int i = head[u]; i; i = pool[i].next) {
                int v = pool[i].v, w = pool[i].w, c = pool[i].c;
                if(w && dis[u] + c < dis[v]) {
                    dis[v] = dis[u] + c;
                    dep[v] = dep[u] + 1;
                    if(vis[v]) continue;
                    vis[v] = 1;
                    que[++tl] = v;
                }
            }
        }
        return dis[t] != INF;
    }
    int dfs(int u, int sum, int cost) {
        if(u == t) {
            ans += sum * cost;
            return sum;
        }
        int flow = 0;
        for(int i = now[u]; i && sum; i = pool[i].next) {
            int v = pool[i].v, w = pool[i].w, c = pool[i].c;
            now[u] = i;
            if(!w || dep[v] != dep[u] + 1 || dis[v] != dis[u] + c || vis[v]) continue;
            int k = dfs(v, min(sum, w), cost + c);
            if(!k) dis[v] = 0;
            else {
                pool[i].w -= k;
                pool[i ^ 1].w += k;
                sum -= k;
                flow += k;
            }
        }
        return flow;
    }
    void calc() {
        int flow = 0;
        while(spfa()) flow += dfs(s, INF, 0);
        cout << flow << ' ' << ans << endl;
    }
}

int main() {

    cin >> n >> m >> s >> t;
    for(int i = 1; i <= m; i++) {
        int u, v, w, c;
        cin >> u >> v >> w >> c;
        addEdge(u, v, w, c);
        addEdge(v, u, 0, -c);
    }

    dinic::calc();

    return 0;
}

原始对偶

但是如果不选择边数最少的增广路，时间复杂度是没有保证的，为 \(O(nmf)\)，非常劣。由于瓶颈是 spfa 不是 dfs，所以我们希望可以将 spfa 用 Dijkstra 代替。但众所周知这个图上面是有负权边的。所以这个地方我们引入势函数 \(h(x)\)，并将每条边的边权置为 \(w(u,v)+h(u)-h(v)\)。发现这样只会将求出来的 \(dis\) 偏移一个常数，所以不会影响最短路。

初始时，令 \(h(x)=d_0(x)\) 即到源点的距离。这个东西需要用一遍 spfa 求出。然后根据三角不等式，所有边权就都变成非负数了，可以直接 Dijkstra。

但是增广可能导致一些边加入，一些边消失，可能影响 \(d(x)\)。但是只需要在每轮增广结束后执行 \(h(x)\gets h(x)+d(x)\)（\(d(x)\) 表示引入势能后的最短路），就能保证加入的反向边也拥有非负边权。我们来证明这一点：注意到新增反向边时当且仅当它位于最短路上，因此满足 \(d(v)=d(u)+(w+h(u)-h(v))\)。此时反向边的边权为 \(w+h'(v)-h'(u)=0\ge 0\)。并且原边的边权显然保持非负。

这样时间复杂度应该是 \(O(n^2f)\) 或 \(O(mf\log n)\)。

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

#include<iostream>
#include<queue>
using namespace std;
const int N = 5010;
const int M = 5e4 + 10;
const int INF = 0x3f3f3f3f;

struct Edge {
    int v, w, c, next;
} pool[2 * M];
int ne = 1, head[N];

void addEdge(int u, int v, int w, int c) {
    pool[++ne] = {v, w, c, head[u]};
    head[u] = ne;
}

int n, m, s, t;

namespace Dinic {
    struct Node {
        int u, d;
        inline bool operator<(const Node &b) const {
            return d > b.d;
        }
    };
    int d[N], now[N], h[N];
    void spfa() {
        static int que[N], hd, tl;
        static bool inq[N];
        for(int i = 1; i <= n; i++) h[i] = INF;
        hd = 1, tl = 0;
        que[++tl] = s;
        inq[s] = 1; h[s] = 0;
        while(hd != (tl + 1) % N) {
            int u = que[hd]; hd = (hd + 1) % N;
            inq[u] = 0;
            for(int e = head[u]; e; e = pool[e].next) {
                int v = pool[e].v, w = pool[e].w, c = pool[e].c;
                if(!w) continue;
                if(h[u] + c < h[v]) {
                    h[v] = h[u] + c;
                    if(!inq[v]) que[(++tl) %= N] = v, inq[v] = 1;
                }
            }
        }
    }
    int dijk() {
        static bool vis[N];
        static priority_queue<Node> que;
        while(!que.empty()) que.pop();
        for(int i = 1; i <= n; i++) now[i] = head[i];
        for(int i = 1; i <= n; i++) vis[i] = 0, d[i] = INF;
        que.push({s, 0}); d[s] = 0;
        while(!que.empty()) {
            int u = que.top().u; que.pop();
            if(vis[u]) continue;
            vis[u] = 1;
            for(int e = head[u]; e; e = pool[e].next) {
                int v = pool[e].v, w = pool[e].w, c = pool[e].c;
                if(!w) continue;
                if(d[u] + (c + h[u] - h[v]) < d[v]) {
                    d[v] = d[u] + (c + h[u] - h[v]);
                    que.push({v, d[v]});
                }
            }
        }
        return d[t] != INF;
    }
    int flow, cost;
    bool vis[N];
    int dfs(int u, int sum) {
        if(u == t) return sum;
        int flow = 0;
        vis[u] = 1;
        for(int e = now[u]; e && sum; e = pool[e].next) {
            now[u] = e;
            int v = pool[e].v, w = pool[e].w, c = pool[e].c;
            if(!w || d[v] != d[u] + (c + h[u] - h[v]) || vis[v]) continue;
            int tmp = dfs(v, min(sum, w));
            if(!tmp) {
                vis[v] = 1;
            } else {
                pool[e].w -= tmp;
                pool[e ^ 1].w += tmp;
                sum -= tmp;
                flow += tmp;
                cost += c * tmp;
            }
        }
        vis[u] = 0;
        return flow;
    }
    void calc() {
        spfa();
        while(dijk()) {
            for(int i = 1; i <= n; i++) vis[i] = 0;
            flow += dfs(s, INF);
            for(int i = 1; i <= n; i++) h[i] += d[i];
        }
        cout << flow << ' ' << cost << '\n';
    }
}

int main() {

    cin >> n >> m >> s >> t;
    for(int i = 1; i <= m; i++) {
        int u, v, w, c; cin >> u >> v >> w >> c;
        addEdge(u, v, w, c);
        addEdge(v, u, 0, -c);
    }

    Dinic::calc();

    return 0;
}