[算法] 最小费用最大流
思想
主要是把无权最大流中找增广路时，每次都增广以费用为计算的最短路，时间复杂度较高，但通常可以计算出 \(n\leq 1e4\) 的数据。
常用定理：每次增广时的最短路长度严格递增，且不存在负环。
Code
几种常见的最小费用最大流的写法：
EK+SPFA
#include <deque>
#include <cstdio>
using namespace std;
const int MAXN = 5e3 + 5;
const int MAXM = 5e4 + 5;
const int INF = 0x7f7f7f7f;
struct Edge { int To, Cap, Cost, Next; };
Edge edge[MAXM << 1];
int head[MAXN], tot = 1;
void Addedge(int u, int v, int w, int c) {
edge[++tot].Next = head[u], edge[tot].To = v, edge[tot].Cap = w, edge[tot].Cost = c, head[u] = tot;
edge[++tot].Next = head[v], edge[tot].To = u, edge[tot].Cap = 0, edge[tot].Cost = -c, head[v] = tot;
}
deque<int> dq;
int dis[MAXN], nxt[MAXN], lim[MAXN];
bool inque[MAXN];
int n, m, s, t, mcmf;
int Min(int x, int y) {
return x < y ? x : y;
}
bool SPFA() {
for (int i = 1; i <= n; i++) dis[i] = INF, lim[i] = INF;
dis[t] = 0, inque[t] = 1, dq.push_back(t);
while (!dq.empty()) {
int u = dq.front(); dq.pop_front(), inque[u] = 0;
for (int i = head[u]; i; i = edge[i].Next) {
int v = edge[i].To;
if (dis[v] > dis[u] - edge[i].Cost && edge[i ^ 1].Cap) {
dis[v] = dis[u] - edge[i].Cost, nxt[v] = i ^ 1, lim[v] = min(lim[u], edge[i ^ 1].Cap);
if (!inque[v]) {
if (!dq.empty() && dis[v] < dis[dq.front()]) dq.push_front(v);
else dq.push_back(v);
}
inque[v] = 1;
}
}
}
return dis[s] != INF;
}
int Update() {
int now = s;
while (now != t) {
int i = nxt[now];
edge[i].Cap -= lim[s];
edge[i ^ 1].Cap += lim[s];
mcmf += lim[s] * edge[i].Cost;
now = edge[i].To;
}
return lim[s];
}
int EK() {
int res = 0;
while (SPFA()) res += Update();
return res;
}
int main() {
scanf("%d %d %d %d", &n, &m, &s, &t);
for (int i = 1, u, v, w, c; i <= m; i++) {
scanf("%d %d %d %d", &u, &v, &w, &c);
Addedge(u, v, w, c);
}
printf("%d ", EK());
printf("%d", mcmf);
return 0;
}
Dinic+SPFA
只是改为了多路增广，但由于最短路的要求过于苛刻，效率与 EK 算法差距不是特别大。
#include <deque>
#include <cstdio>
using namespace std;
const int MAXN = 5e3 + 5;
const int MAXM = 5e4 + 5;
const int INF = 0x7f7f7f7f;
struct Edge { int To, Cap, Cost, Next; };
Edge edge[MAXM << 1];
int head[MAXN], tot = 1;
void Addedge(int u, int v, int w, int c) {
edge[++tot].Next = head[u], edge[tot].To = v, edge[tot].Cap = w, edge[tot].Cost = c, head[u] = tot;
edge[++tot].Next = head[v], edge[tot].To = u, edge[tot].Cap = 0, edge[tot].Cost = -c, head[v] = tot;
}
deque<int> dq;
int dis[MAXN];
bool inque[MAXN], vis[MAXN];
int n, m, s, t, mcmf;
int Min(int x, int y) {
return x < y ? x : y;
}
bool SPFA() {
for (int i = 1; i <= n; i++) dis[i] = INF;
dis[t] = 0, inque[t] = 1, dq.push_back(t);
while (!dq.empty()) {
int u = dq.front(); dq.pop_front(), inque[u] = 0;
for (int i = head[u]; i; i = edge[i].Next) {
int v = edge[i].To;
if (dis[v] > dis[u] - edge[i].Cost && edge[i ^ 1].Cap) {
dis[v] = dis[u] - edge[i].Cost;
if (!inque[v]) {
if (!dq.empty() && dis[v] < dis[dq.front()]) dq.push_front(v);
else dq.push_back(v);
}
inque[v] = 1;
}
}
}
return dis[s] != INF;
}
int dfs(int u, int flow) {
if (!flow || u == t) return flow;
int rest = flow;
vis[u] = 1;
for (int i = head[u]; i && rest; i = edge[i].Next) {
int v = edge[i].To;
if (!vis[v] && dis[v] == dis[u] - edge[i].Cost && edge[i].Cap) {
int del = dfs(v, Min(edge[i].Cap, rest));
edge[i].Cap -= del, edge[i ^ 1].Cap += del, rest -= del;
mcmf += del * edge[i].Cost;
}
}
vis[u] = 0;
return flow - rest;
}
int Dinic() {
int res = 0, flow;
while (SPFA()) while ((flow = dfs(s, INF))) res += flow;
return res;
}
int main() {
scanf("%d %d %d %d", &n, &m, &s, &t);
for (int i = 1, u, v, w, c; i <= m; i++) {
scanf("%d %d %d %d", &u, &v, &w, &c);
Addedge(u, v, w, c);
}
printf("%d ", Dinic());
printf("%d", mcmf);
return 0;
}
EK+Dijkstra
将 SPFA 改为带势 Dijkstra，将 \(w[u][v]\) 转化为 \(h[u]-h[v]+w[u][v]\)，那么 \(s\to t\) 的中间节点的势函数都被消掉了。
但需要保证 \(h[u]-h[v]+w[u][v]\geq0\) ，对于每次求最短路时 \(dis\) 不为极大值的点，势函数加上最短路。即 \(h[u]+=dis[u](dis[u]\neq \infty)\)。
证明该构造方案合法：
\(h[u]-h[v]+w[u][v]\geq0\) 移项变为 \(h[u]+w[u][v]\geq h[v]\)。
- \(w[u][v]>0\) 时，显然成立。
- \(w[u][v]<0\) 。若当前反边是新增的，那么存在一次正边被增广。即 \(dis'[v]-w[u][v]=dis'[u]\to dis'[v]=dis'[u]+w[u][v]\)，由于正边是在情况 \(1\) 中成立的，那么此时的 \(dis'[v]<dis'[u]\)，而累加到 \(h\) 中刚好将 \(w[u][v]\) 抵消了，则也成立，简单来说，就是正边的增广将反边的增广相抵消了。否则，\(h=\Delta dis+k\) ，而当前边的方向并没有改变，则 \(\Delta dis\) 并不会发生改变，其中 \(k\) 是由最短路中的其他边所构成的，\(u,v\) 的 \(k\) 相等，并不会影响正确性。
若图中存在负边，那么先用 SPFA 求出第一次的势函数，再用 Dijkstra 修改势函数。
#include <deque>
#include <cstdio>
using namespace std;
const int MAXN = 5e3 + 5;
const int MAXM = 5e4 + 5;
const int INF = 0x7f7f7f7f;
struct Edge { int To, Cap, Cost, Next; };
Edge edge[MAXM << 1];
int head[MAXN], tot = 1;
void Addedge(int u, int v, int w, int c) {
edge[++tot].Next = head[u], edge[tot].To = v, edge[tot].Cap = w, edge[tot].Cost = c, head[u] = tot;
edge[++tot].Next = head[v], edge[tot].To = u, edge[tot].Cap = 0, edge[tot].Cost = -c, head[v] = tot;
}
deque<int> dq;
int dis[MAXN], nxt[MAXN], lim[MAXN];
bool inque[MAXN];
int n, m, s, t, mcmf;
int Min(int x, int y) {
return x < y ? x : y;
}
bool SPFA() {
for (int i = 1; i <= n; i++) dis[i] = INF, lim[i] = INF;
dis[t] = 0, inque[t] = 1, dq.push_back(t);
while (!dq.empty()) {
int u = dq.front(); dq.pop_front(), inque[u] = 0;
for (int i = head[u]; i; i = edge[i].Next) {
int v = edge[i].To;
if (dis[v] > dis[u] - edge[i].Cost && edge[i ^ 1].Cap) {
dis[v] = dis[u] - edge[i].Cost, nxt[v] = i ^ 1, lim[v] = min(lim[u], edge[i ^ 1].Cap);
if (!inque[v]) {
if (!dq.empty() && dis[v] < dis[dq.front()]) dq.push_front(v);
else dq.push_back(v);
}
inque[v] = 1;
}
}
}
return dis[s] != INF;
}
int Update() {
int now = s;
while (now != t) {
int i = nxt[now];
edge[i].Cap -= lim[s];
edge[i ^ 1].Cap += lim[s];
mcmf += lim[s] * edge[i].Cost;
now = edge[i].To;
}
return lim[s];
}
int EK() {
int res = 0;
while (SPFA()) res += Update();
return res;
}
int main() {
scanf("%d %d %d %d", &n, &m, &s, &t);
for (int i = 1, u, v, w, c; i <= m; i++) {
scanf("%d %d %d %d", &u, &v, &w, &c);
Addedge(u, v, w, c);
}
printf("%d ", EK());
printf("%d", mcmf);
return 0;
}
例：[SDOI2017]新生舞会
题意
有 \(n\) 个点，第 \(i\) 个点与第 \(j\) 个点匹配会产生 \(a_{i,j}\) 和 \(b_{i,j}\)，每个点只能匹配一次。
求 \(\frac{\sum a}{\sum b}\) 的最大值。
思路
二分图带权匹配，再用 \(0/1\) 分数规划的套路，使用费用流求解即可
Code
#include <map>
#include <queue>
#include <cstdio>
#include <cstring>
using namespace std;
const int MAXN = 2e2 + 5;
const int MAXM = 2e4 + 5;
const int INF = 0x3f3f3f3f;
#define Eps 1e-8
struct Edge { int To, Cap, Next; long double Cost; };
Edge edge[MAXM << 1];
int head[MAXN], tot = 1;
void Addedge(int u, int v, int w, long double c) {
edge[++tot].Next = head[u], edge[tot].To = v, edge[tot].Cap = w, edge[tot].Cost = c, head[u] = tot;
edge[++tot].Next = head[v], edge[tot].To = u, edge[tot].Cap = 0, edge[tot].Cost = -c, head[v] = tot;
}
priority_queue<pair<int, int>, vector<pair<int, int> >, greater<pair<int, int> > > pq;
long double dis[MAXN], h[MAXN], mcmf;
int dep[MAXN], pre[MAXN], lim[MAXN];
bool vis[MAXN];
int a[105][105], b[105][105];
int n, m, s, t;
int Min(int x, int y) {
return x < y ? x : y;
}
bool Dijkstra() {
for (int i = s; i <= t; i++) vis[i] = 0, pre[i] = 0, lim[i] = dis[i] = INF;
dis[s] = 0, pq.push(make_pair(0, s));
while (!pq.empty()) {
int u = pq.top().second; pq.pop();
if (vis[u]) continue; vis[u] = 1;
for (int i = head[u]; i; i = edge[i].Next) {
int v = edge[i].To;
if (dis[v] > dis[u] + edge[i].Cost + h[u] - h[v] && edge[i].Cap) {
pre[v] = i, lim[v] = min(lim[u], edge[i].Cap);
dis[v] = dis[u] + edge[i].Cost + h[u] - h[v];
pq.push(make_pair(dis[v], v));
}
}
}
for (int i = s; i <= t; i++) if (dis[i] < INF) h[i] += dis[i];
return dis[t] != INF;
}
int Update() {
int now = t;
while (now != s) {
int i = pre[now];
edge[i].Cap -= lim[t];
edge[i ^ 1].Cap += lim[t];
mcmf += lim[t] * edge[i].Cost;
now = edge[i ^ 1].To;
}
return lim[t];
}
int EK() {
int res = 0;
while (Dijkstra()) res += Update();
return res;
}
bool Check(long double mid) {
for (int i = s; i <= t; i++) head[i] = 0, h[i] = 0; tot = 1; mcmf = 0;
t = 2 * n + 1;
for (int i = 1; i <= n; i++) Addedge(s, i, 1, 0);
for (int i = 1; i <= n; i++) Addedge(i + n, t, 1, 0);
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
Addedge(i, j + n, 1, mid * b[i][j] - a[i][j]);
}
}
EK();
return mcmf <= 0;
}
int main() {
scanf("%d", &n);
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
scanf("%d", &a[i][j]);
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
scanf("%d", &b[i][j]);
long double l = 0, r = 1e4;
while (r - l >= Eps) {
long double mid = (l + r) / 2;
if (Check(mid)) l = mid;
else r = mid;
}
printf("%.6Lf", l);
return 0;
}
