[算法] 最小费用最大流

思想

主要是把无权最大流中找增广路时，每次都增广以费用为计算的最短路，时间复杂度较高，但通常可以计算出 \(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;
}