费用流之 Cost Scaling Algorithm

__EternalLife__ · 2026-10-14 20:35:21 · 算法·理论

前言

起因是模拟赛搬了 AT_agc018_c，大致题意如下：

有 x+y+z 件物品，每件物品有三种价值 a_i,b_i,c_i，你要选择 x 件物品价值为 a，y 件物品价值为 b，z 件物品价值为 c，求价值之和的最大值。

第一眼转化为费用流问题，然后使用了高贵的 SSP 算法，于是荣获 TLE。

信奉网络流的我坚定地认为费用流可以创过这个题，于是通过赛后拷问 AI 发现了这个算法，即成本缩放算法（Cost Scaling Algorithm），然后成功创过了此题，遂有此文。

后面都用 CS 算法代称。

介绍

CS 算法一般用于求解最小费用循环流和最小费用固定流问题，其核心思路是先允许一定费用误差，然后逐步缩小误差。

一些定义与引入

超额流量

对于每个点 v，定义超额流量 ex_v 表示流入 v 的流量减去流出 v 的流量，则

如果最终的 ex_v 都等于 0，则说明流合法。

如果要求解的是最小费用固定流问题，则需初始化 ex_s=F,ex_t=-F。

势能与约化费用

由于 CS 需要通过缩小误差来得到答案，那么调整费用肯定是不可避免的，为了避免调整对答案带来的影响，我们引入势能 p_v，然后定义约化费用为

这样，在后续调整中，势能会改变单条边和路径的约化费用，但同起点、同终点的路径之间，约化费用差等于原费用差，且环费用完全不变。

\varepsilon- 最优

定义如果残量网络中的每条边都满足

那么称当前流是 \varepsilon- 最优，最终达到 \varepsilon<\frac{1}{n} 时，我们便找到了最优方案。 :::info[证明] 假设最终流不是最优的。那么一定存在一个负费用残量环。考虑任意一个简单残留环。由于当前流是 \varepsilon- 最优的，故每条边满足

所以整个环的约化费用满足

这样一次可以推很多流，因为只要 u 有多余流量，边还有残量且边的约化费用为负，那么沿它推任何不超过上述限制的流量，都会让费用下降。因此不需要像最短路增广那样每次只发送一条路径上的瓶颈流，这也是 CS 算法高效之处。

Relabel

假设某个点 u 满足 ex_u>0，且它的所有残量出边都满足 \bar{c}(u,v)\geq 0，则此时只能改变 u 的势能，使得有边可以推进，即

以上两个过程称为一次 refine，具体实现如下：

ll get(int u,int i){
    return e[i].cost+p[u]-p[e[i].to];
}
void push(int u,int i,ll f){
    e[i].cap-=f; e[i^1].cap+=f;
    ex[u]-=f; ex[e[i].to]+=f; 
}
void refine(ll eps){
    for(int u=1;u<=n;u++){
        for(int i=head[u];i;i=e[i].next){
            if(e[i].cap>0&&get(u,i)<0) push(u,i,e[i].cap);
        }
    }
    fill(vis+1,vis+n+1,0);
    queue<int> q;
    for(int u=1;u<=n;u++){
        if(ex[u]>0) q.push(u),vis[u]=1;
    }
    memcpy(cur,head,sizeof(int)*(n+1));
    while(!q.empty()){
        int u=q.front(); q.pop();
        vis[u]=0;
        while(ex[u]>0){
            bool pushed=0;
            while(cur[u]){
                int i=cur[u];
                if(e[i].cap>0&&get(u,i)<0){//有推进边，则Push
                    int v=e[i].to;
                    ll od=ex[v];
                    ll f=min(ex[u],e[i].cap);
                    push(u,i,f);
                    pushed=1;
                    if(od<=0&&ex[v]>0&&!vis[v]) q.push(v),vis[v]=1;
                    if(!e[i].cap) cur[u]=e[i].next;
                    if(!ex[u]) break;
                }
                else cur[u]=e[i].next;
            }
            if(!ex[u]) break;
            if(!pushed){//若无可推进边，则Relabel
                ll mi=LLONG_MAX;
                for(int i=head[u];i;i=e[i].next){
                    if(e[i].cap>0) mi=min(mi,get(u,i));
                }
                p[u]-=mi+eps; cur[u]=head[u];
            }
        }
    }
}

很好理解吧。

减小 \varepsilon

当一个阶段经过 Push 和 Relabel 操作后，满足任意 ex_u=0，则缩小 \varepsilon，一般都是减半处理。当然你减多少都可以，但不一定会使得代码效率更高，因为有可能要做更多 Push 和 Relabel。

时间复杂度

时间复杂度是 O(n^2m\log (nC))。好像还能用动态树优化到 O(nm\log(\frac{n^2}{m})\log(nC))，不过我不会。

解决最小费用最大流问题

由于 CS 是用于求解最小费用循环流问题的，以模板题为例，要想求解最小费用最大流，有两个做法：

1. 先用最大流算法求出最大流，然后跑 CS。 :::success[Code from GPT 1.19s]

   
   #include <bits/stdc++.h>
   using namespace std;

using int64 = long long; using i128 = __int128_t;

// 第一步：Dinic，只求最大流流量 struct Dinic { struct Edge { int to, rev; int64 cap; };

int n;
vector<vector<Edge>> g;
vector<int> level, cur;

Dinic(int n) : n(n), g(n), level(n), cur(n) {}

void addEdge(int u, int v, int64 cap) {
    Edge a{v, (int)g[v].size(), cap};
    Edge b{u, (int)g[u].size(), 0};
    g[u].push_back(a);
    g[v].push_back(b);
}

bool bfs(int s, int t) {
    fill(level.begin(), level.end(), -1);
    queue<int> q;

    level[s] = 0;
    q.push(s);

    while (!q.empty()) {
        int u = q.front();
        q.pop();

        for (const Edge &e : g[u]) {
            if (e.cap > 0 && level[e.to] == -1) {
                level[e.to] = level[u] + 1;
                q.push(e.to);
            }
        }
    }

    return level[t] != -1;
}

int64 dfs(int u, int t, int64 f) {
    if (u == t) return f;

    for (int &i = cur[u]; i < (int)g[u].size(); ++i) {
        Edge &e = g[u][i];

        if (e.cap == 0 || level[e.to] != level[u] + 1) {
            continue;
        }

        int64 pushed = dfs(e.to, t, min(f, e.cap));
        if (pushed == 0) continue;

        e.cap -= pushed;
        g[e.to][e.rev].cap += pushed;
        return pushed;
    }

    return 0;
}

int64 maxFlow(int s, int t) {
    int64 flow = 0;
    const int64 INF = numeric_limits<int64>::max() / 4;

    while (bfs(s, t)) {
        fill(cur.begin(), cur.end(), 0);
        while (int64 pushed = dfs(s, t, INF)) {
            flow += pushed;
        }
    }

    return flow;
}

};

// 第二步：Cost Scaling，求恰好 F 单位流的最小费用 struct CostScaling { struct Edge { int to, rev; int64 cap, cost; };

int n;
vector<vector<Edge>> g;

// pot[v]：点势能
vector<int64> pot;

// excess[v]：当前点的超额流量
vector<int64> excess;

// 缩放后的总费用
i128 totalCost = 0;

CostScaling(int n) : n(n), g(n), pot(n, 0), excess(n, 0) {}

void addEdge(int u, int v, int64 cap, int64 cost) {
    Edge forward{v, (int)g[v].size(), cap, cost};
    Edge backward{u, (int)g[u].size(), 0, -cost};

    g[u].push_back(forward);
    g[v].push_back(backward);
}

int64 reducedCost(int u, const Edge &e) const {
    return e.cost + pot[u] - pot[e.to];
}

// 沿 u 的第 i 条残量边推送 f 单位流量
void push(int u, int i, int64 f) {
    Edge &e = g[u][i];
    int v = e.to;
    int rev = e.rev;

    e.cap -= f;
    g[v][rev].cap += f;

    excess[u] -= f;
    excess[v] += f;

    totalCost += (i128)f * e.cost;
}

// 求恰好 F 单位流的最小费用
bool sendFixedFlow(int s, int t, int64 F) {
    int64 maxCost = 0;
    for (const auto &adj : g) {
        for (const Edge &e : adj) {
            maxCost = max(maxCost, llabs(e.cost));
        }
    }

    // 初始 epsilon 取不小于最大费用的最小二次幂
    int64 eps = 1;
    while (eps < maxCost) {
        eps <<= 1;
    }

    // excess = 外部供给 + 流入 - 流出
    excess[s] += F;
    excess[t] -= F;

    while (eps > 0) {
        // 饱和所有负约化费用残量边
        for (int u = 0; u < n; ++u) {
            for (int i = 0; i < (int)g[u].size(); ++i) {
                if (g[u][i].cap > 0 &&
                    reducedCost(u, g[u][i]) < 0) {
                    push(u, i, g[u][i].cap);
                }
            }
        }

        // 把所有有正超额的点放入队列
        queue<int> q;
        vector<char> inQueue(n, false);

        for (int v = 0; v < n; ++v) {
            if (excess[v] > 0) {
                q.push(v);
                inQueue[v] = true;
            }
        }

        // 不断消除正超额：Push 或 Relabel
        while (!q.empty()) {
            int u = q.front();
            q.pop();
            inQueue[u] = false;

            int cur = 0;

            while (excess[u] > 0) {
                bool pushedAny = false;

                // 寻找约化费用为负的可推进边
                while (cur < (int)g[u].size()) {
                    Edge &e = g[u][cur];

                    if (e.cap > 0 && reducedCost(u, e) < 0) {
                        int v = e.to;
                        int64 amount = min(excess[u], e.cap);
                        int64 oldExcessV = excess[v];

                        push(u, cur, amount);
                        pushedAny = true;

                        // v 从非正超额变成正超额时，加入队列
                        if (oldExcessV <= 0 &&
                            excess[v] > 0 &&
                            !inQueue[v]) {
                            q.push(v);
                            inQueue[v] = true;
                        }

                        if (g[u][cur].cap == 0) {
                            ++cur;
                        }

                        if (excess[u] == 0) {
                            break;
                        }
                    } else {
                        ++cur;
                    }
                }

                if (excess[u] == 0) {
                    break;
                }

                // 没有可推进边：重标记
                if (!pushedAny) {
                    int64 minRC = numeric_limits<int64>::max();

                    for (const Edge &e : g[u]) {
                        if (e.cap > 0) {
                            minRC = min(minRC, reducedCost(u, e));
                        }
                    }

                    // F 是最大流阶段得到的可行流量，正常情况下不会触发
                    if (minRC == numeric_limits<int64>::max()) {
                        return false;
                    }

                    // 将最便宜的残量出边调整到约化费用 -eps
                    pot[u] -= minRC + eps;
                    cur = 0;
                }
            }
        }

        eps >>= 1;
    }

    for (int64 x : excess) {
        if (x != 0) {
            return false;
        }
    }

    return true;
}

};

void printInt128(i128 x) { if (x == 0) { cout << '0'; return; }

if (x < 0) {
    cout << '-';
    x = -x;
}

string s;
while (x > 0) {
    s.push_back(char('0' + x % 10));
    x /= 10;
}

reverse(s.begin(), s.end());
cout << s;

}

int main() { ios::sync_with_stdio(false); cin.tie(nullptr);

int n, m, s, t;
cin >> n >> m >> s >> t;
--s;
--t;

struct InputEdge {
    int u, v;
    int64 cap, cost;
};

vector<InputEdge> edges;
edges.reserve(m);

Dinic dinic(n);

for (int i = 0; i < m; ++i) {
    int u, v;
    int64 cap, cost;
    cin >> u >> v >> cap >> cost;
    --u;
    --v;

    edges.push_back({u, v, cap, cost});
    dinic.addEdge(u, v, cap);
}

// 求最大流量 F
int64 maxFlow = dinic.maxFlow(s, t);

// 把费用放大 n+1 倍，使用整数 epsilon 缩放
const int64 SCALE = n + 1;

CostScaling costScaling(n);
for (const auto &e : edges) {
    costScaling.addEdge(
        e.u,
        e.v,
        e.cap,
        e.cost * SCALE
    );
}

// 求恰好 maxFlow 单位流的最小费用
if (!costScaling.sendFixedFlow(s, t, maxFlow)) {
    // 理论上不会发生：maxFlow 已经保证这个流量可行
    return 0;
}

cout << maxFlow << ' ';
printInt128(costScaling.totalCost / SCALE);
cout << '\n';

return 0;

}

:::
2. 连一条从 $t$ 向 $s$ 的边，$cap(t,s)=\sum_e cap_e$，$cost(t,s)=-(\sum_e cap_e|cost_e|+1))$。
:::success[我的大常实现 2.92s]

include<bits/stdc++.h>

define endl '\n'

using namespace std; using ll=long long; const int maxn=1e5+10; int n,m,s,t; struct node{ int to,next; ll cap,cost; }e[maxn2]; int head[maxn],cur[maxn],tot=1; ll p[maxn],ex[maxn]; void adde(int u,int v,ll cap,ll cost){ e[++tot]={v,head[u],cap,cost}; head[u]=tot; } void add(int u,int v,ll cap,ll cost){ adde(u,v,cap,cost); adde(v,u,0,-cost); } ll get(int u,int i){ return e[i].cost+p[u]-p[e[i].to]; } void push(int u,int i,ll f){ e[i].cap-=f; e[i^1].cap+=f; ex[u]-=f; ex[e[i].to]+=f; } bool vis[maxn]; void refine(ll eps){ for(int u=1;u<=n;u++){ for(int i=head[u];i;i=e[i].next){ if(e[i].cap>0&&get(u,i)<0) push(u,i,e[i].cap); } } fill(vis+1,vis+n+1,0); queue<int> q; for(int u=1;u<=n;u++){ if(ex[u]>0) q.push(u),vis[u]=1; } memcpy(cur,head,sizeof(int)(n+1)); while(!q.empty()){ int u=q.front(); q.pop(); vis[u]=0; while(ex[u]>0){ bool pushed=0; while(cur[u]){ int i=cur[u]; if(e[i].cap>0&&get(u,i)<0){ int v=e[i].to; ll od=ex[v]; ll f=min(ex[u],e[i].cap); push(u,i,f); pushed=1; if(od<=0&&ex[v]>0&&!vis[v]) q.push(v),vis[v]=1; if(!e[i].cap) cur[u]=e[i].next; if(!ex[u]) break; } else cur[u]=e[i].next; } if(!ex[u]) break; if(!pushed){ ll mi=LLONG_MAX; for(int i=head[u];i;i=e[i].next){ if(e[i].cap>0) mi=min(mi,get(u,i)); } p[u]-=mi+eps; cur[u]=head[u]; } } } } void costScaling(){ ll mx=1; for(int i=2;i<=tot;i++) mx=max(mx,llabs(e[i].cost)); ll eps=1; while(eps<mx) eps<<=1; while(eps>0) refine(eps),eps>>=1; } int u[maxn],v[maxn]; ll cap[maxn],cost[maxn]; signed main(){ ios::sync_with_stdio(0); cin.tie(0); cout.tie(0); cin>>n>>m>>s>>t; ll sum1=0,sum2=0; for(int i=1;i<=m;i++) cin>>u[i]>>v[i]>>cap[i]>>cost[i],sum1+=cap[i],sum2+=cap[i]abs(cost[i]); vector<int> id; for(int i=1;i<=m;i++) id.push_back(tot+1),add(u[i],v[i],cap[i],cost[i](n+1)); int ed=tot+1; add(t,s,sum1+1,-(sum2+1)(n+1)); costScaling(); cout<<sum1+1-e[ed].cap<<" "; ll c=0; for(int i=1;i<=m;i++){ ll used=cap[i]-e[id[i-1]].cap; c+=usedcost[i]; } cout<<c; return 0; } ::: :::success[GPT 的小常实现 1.49s]

#include <bits/stdc++.h>
using namespace std;

using int64 = long long;
using i128 = __int128_t;

struct CostScalingCirculation {
    struct Edge {
        int to;
        int rev;
        int64 cap;   // 残量容量
        int64 cost;  // 缩放后的费用
    };

    int n;
    vector<vector<Edge>> g;

    // pot[v]：点势能
    vector<int64> pot;

    // excess[v]：点 v 当前的超额流量
    vector<int64> excess;

    // 当前弧
    vector<int> cur;

    CostScalingCirculation(int n)
        : n(n), g(n), pot(n, 0), excess(n, 0), cur(n, 0) {}

    // 加入正向边和反向残量边，返回正向边在 g[u] 中的下标
    int addEdge(int u, int v, int64 cap, int64 cost) {
        int index = (int)g[u].size();

        Edge forward{v, (int)g[v].size(), cap, cost};
        Edge backward{u, index, 0, -cost};

        g[u].push_back(forward);
        g[v].push_back(backward);

        return index;
    }

    // 约化费用：c(u,v) + pot[u] - pot[v]
    int64 reducedCost(int u, const Edge &e) const {
        return e.cost + pot[u] - pot[e.to];
    }

    // 沿 u 的第 i 条边推出 f 单位流量
    void push(int u, int i, int64 f) {
        Edge &e = g[u][i];
        int v = e.to;
        int rev = e.rev;

        e.cap -= f;
        g[v][rev].cap += f;

        excess[u] -= f;
        excess[v] += f;
    }

    // 执行一个 epsilon 阶段
    void refine(int64 epsilon) {
        /*
         * 先饱和所有约化费用小于 0 的残量边。
         * 这会产生伪流和 excess，之后再通过 Push/Relabel 消除。
         */
        for (int u = 0; u < n; ++u) {
            for (int i = 0; i < (int)g[u].size(); ++i) {
                if (g[u][i].cap > 0 &&
                    reducedCost(u, g[u][i]) < 0) {
                    int64 f = g[u][i].cap;
                    push(u, i, f);
                }
            }
        }

        queue<int> q;
        vector<char> inQueue(n, false);

        for (int v = 0; v < n; ++v) {
            if (excess[v] > 0) {
                q.push(v);
                inQueue[v] = true;
            }
        }

        fill(cur.begin(), cur.end(), 0);

        /*
         * 不断处理正 excess 点：
         * 有负约化费用边就 Push；
         * 否则 Relabel，让最便宜的残量出边变为可推进边。
         */
        while (!q.empty()) {
            int u = q.front();
            q.pop();
            inQueue[u] = false;

            while (excess[u] > 0) {
                bool pushed = false;

                while (cur[u] < (int)g[u].size()) {
                    int i = cur[u];
                    Edge &e = g[u][i];

                    if (e.cap > 0 && reducedCost(u, e) < 0) {
                        int v = e.to;
                        int64 oldExcess = excess[v];
                        int64 f = min(excess[u], e.cap);

                        push(u, i, f);
                        pushed = true;

                        // v 从非正 excess 变成正 excess 时，加入队列
                        if (oldExcess <= 0 &&
                            excess[v] > 0 &&
                            !inQueue[v]) {
                            q.push(v);
                            inQueue[v] = true;
                        }

                        if (g[u][i].cap == 0) {
                            ++cur[u];
                        }

                        if (excess[u] == 0) {
                            break;
                        }
                    } else {
                        ++cur[u];
                    }
                }

                if (excess[u] == 0) {
                    break;
                }

                if (!pushed) {
                    int64 minReducedCost = LLONG_MAX;

                    for (const Edge &e : g[u]) {
                        if (e.cap > 0) {
                            minReducedCost =
                                min(minReducedCost, reducedCost(u, e));
                        }
                    }

                    // 本题转化出的循环流总是可行，此断言应成立
                    assert(minReducedCost != LLONG_MAX);

                    /*
                     * 将最便宜的残量出边约化费用调整到 -epsilon。
                     * 只修改势能，不改变流量。
                     */
                    pot[u] -= minReducedCost + epsilon;
                    cur[u] = 0;
                }
            }
        }
    }

    // 求最小费用循环流
    void minCostCirculation() {
        int64 maxAbsCost = 1;

        for (const auto &adj : g) {
            for (const Edge &e : adj) {
                maxAbsCost = max(maxAbsCost, llabs(e.cost));
            }
        }

        // 初始 epsilon 取不小于最大边费用绝对值的最小二次幂
        int64 epsilon = 1;
        while (epsilon < maxAbsCost) {
            epsilon <<= 1;
        }

        while (epsilon > 0) {
            refine(epsilon);
            epsilon >>= 1;
        }

        for (int64 x : excess) {
            assert(x == 0);
        }
    }
};

void printInt128(i128 x) {
    if (x == 0) {
        cout << '0';
        return;
    }

    if (x < 0) {
        cout << '-';
        x = -x;
    }

    string s;
    while (x > 0) {
        s.push_back(char('0' + x % 10));
        x /= 10;
    }

    reverse(s.begin(), s.end());
    cout << s;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int n, m, s, t;
    cin >> n >> m >> s >> t;
    --s;
    --t;

    struct InputEdge {
        int u, v;
        int64 cap, cost;
    };

    vector<InputEdge> edges;
    edges.reserve(m);

    int64 sumCapacity = 0;
    i128 totalCostBound128 = 0;

    for (int i = 0; i < m; ++i) {
        int u, v;
        int64 cap, cost;
        cin >> u >> v >> cap >> cost;
        --u;
        --v;

        edges.push_back({u, v, cap, cost});
        sumCapacity += cap;
        totalCostBound128 += (i128)cap * abs(cost);
    }
    int64 totalCostBound = (int64)totalCostBound128;
    int64 penalty = totalCostBound + 1;
    const int64 SCALE = n + 1;

    CostScalingCirculation cs(n);

    // 记录每条原图边在残量网络中的位置，最后据此统计费用
    vector<pair<int, int>> edgeRef;
    edgeRef.reserve(m);

    for (const auto &e : edges) {
        int index = cs.addEdge(
            e.u, e.v, e.cap, e.cost * SCALE
        );
        edgeRef.push_back({e.u, index});
    }

    /*
     * 添加人工边 t -> s。
     * 其容量大于任何可能的原图流量。
     */
    int64 artificialCapacity = sumCapacity + 1;
    int artificialEdgeIndex = cs.addEdge(
        t, s, artificialCapacity, -penalty * SCALE
    );

    // 纯 Cost Scaling 求最小费用循环流
    cs.minCostCirculation();

    // 人工边上的流量就是原图的最大流量
    int64 maxFlow =
        artificialCapacity - cs.g[t][artificialEdgeIndex].cap;

    // 根据每条原图正向边的残量，恢复原图流量并统计原始费用
    i128 minCost = 0;

    for (int i = 0; i < m; ++i) {
        int u = edgeRef[i].first;
        int index = edgeRef[i].second;

        int64 usedFlow =
            edges[i].cap - cs.g[u][index].cap;

        minCost += (i128)usedFlow * edges[i].cost;
    }

    cout << maxFlow << ' ';
    printInt128(minCost);
    cout << '\n';

    return 0;
}

::: 由于此算法优良的性质，对于有负圈的费用流直接做即可，不需要做什么复杂的处理，而且跑得要比 SSP 快得多。

对比