最小费用最大流

1.问题

给定一张图，有 $n$ 个节点，给定源点和汇点，有 $m$ 条边，给定起点 $u$ 与终点 $v$ 并给出其容量 $w$ 与流过单位流量需耗费的费用 $x$ ，求最小费用。

2.解决方案

将每条边的单位流量看作其边权，将原本的 $Dinic$ 的 $bfs$ 分层改为$SPFA/Dij$$kstra$ 求最短路，每次增广最短路，求出的显然是最大流，又因为最大流算法可以回流，故求出的也一定是最小费用。
但如何保证在此过程中不会出现负环，导致我们的最短路算法炸掉呢?

3.证明

首先我们证明这样一个问题：

增广路的长度一定非严格递增。

尝试用反证法：

假设出现了在第$i$次时这种情况，那么一定是在第 $i-1$ 次增广后出现了更短的最短路。我们不妨假设新出现的边是 $v$ 到 $u$ 。

如下图，$u$ 到 $v$ 有一条费用为 $f$ 的边，新出现的 $v$ 到 $u$ 的边费用为$-f$，$s$ 到 $u$ 的路径费用为 $a$，$u$ 到 $t$ 的路径费用为 $b$ ；$s$到$v$ 的路径费用为 $c$，$v$ 到 $t$ 的路径费用为 $d$。

由于 $i-1$ 次增广走了 $s->u->v->t$。
显然有 $a+f<=c$，$f+d<=b$。
所以 $b+c-f>=a+d+f$。
也就是说第$i$次增广的最短路长度不可能比第$i-1$次的短。
故可以证明出增广路长度非严格递增，同时也证明了在增广过程中不可能出现负环。

4.代码

下面给出$SPFA+Dinic$的模板：

#include<bits/stdc++.h>
using namespace std;
long long T,s,t,n,m,k,ans,ans1;
long long d[10005],cur[10005],v[2000005],ne[2000005],h[10005],val[2000005],val1[2000005],cnt=1,p[10005];
bool f[10005];
void add(long long x,long long y,long long z,long long a){//加双向边
	v[++cnt]=y;
	ne[cnt]=h[x];
	val[cnt]=z;
	val1[cnt]=a;
	h[x]=cnt;
	swap(x,y),z=0;
	v[++cnt]=y;
	ne[cnt]=h[x];
	val[cnt]=0;
	val1[cnt]=-a;//费用为负
	h[x]=cnt;
}
bool bfs(){//SPFA
	memset(d,0x3f,sizeof d);//初始化
	memset(f,0,sizeof f);
	for(int i=1;i<=n;i++)cur[i]=h[i];
	queue<int>q;
	d[s]=0;
	f[s]=1;
	q.push(s);
	while(!q.empty()){
		long long u=q.front();
		q.pop();
		f[u]=0;
		for(long long i=h[u];i;i=ne[i]){
			if(val[i]==0)continue;
			long long vv=v[i],w=val1[i];
			if(d[vv]>d[u]+w){
				d[vv]=d[u]+w;
				if(!f[vv]){
					f[vv]=1;
					q.push(vv);
				}
			}
		}
	}
	memset(f,0,sizeof f);
	return d[t]!=0x3f3f3f3f3f3f3f3f;
}
long long dfs(long long x,long long y){
	if(x==t||y==0){
		return y;
	}
	f[x]=1;
	long long rp=0;
	for(long long i=cur[x];i;i=ne[i]){
		if(f[v[i]]==0&&d[v[i]]==d[x]+val1[i]&&val[i]){
			cur[x]=i;//弧优化
			long long kk=dfs(v[i],min(y-rp,val[i]));
			rp+=kk;
			val[i]-=kk;
			val[i^1]+=kk;
			ans1+=kk*val1[i];//统计费用
			if(y==rp)break;
		}
	}
	f[x]=0;
	return rp;
}
int main(){
	scanf("%lld%lld%lld%lld",&n,&m,&s,&t);
	for(long long i=1,u,v,w,x;i<=m;i++){
		scanf("%lld%lld%lld%lld",&u,&v,&w,&x);
		add(u,v,w,x);//加边
	}
	while(bfs())ans+=dfs(s,0x3f3f3f3f3f3f3f3f3f);
	printf("%lld %lld",ans,ans1);
}

5.关于Dijkstra

因为 SPFA 已死，故我们要用Dijkstra来求解最短路。

$Dijkstra$ 由于不能跑有负边权的图，
故引入一个新的东西:势能函数 $h_i$。
它用来保证边权也就是 $w_{i,j}+h_i-h_j$ 不为负，
在一条路径上，h抵消后只会剩下 $h_s-h_v$，$v$ 为路径终点。
在求解完最短路后只用减掉 $h_s-h_v$ 即可.
那么如何设置 $h$ 呢?
我们可以将h设置为上一次求解的 $s$ 到每一个点的最短路 $dist$ ，
又通过上面证过的结论：增广路长度非严格递增，
我们可以得知 $w_{i,j}+h_i-h_j>=0$。
总而言之：
我们只在开始时用 $SPFA$ 求解一遍 $dist$ 。
之后便用 $dist$ 更新 $h$ 即可。
代码：

#include<bits/stdc++.h>
using namespace std;
long long T,s,t,n,m,k,ans,ans1;
long long d[10005],sn[10005],cur[10005],v[2000005],ne[2000005],h[10005],val[2000005],val1[2000005],cnt=1;
bool f[10005];
void add(long long x,long long y,long long z,long long a){
	v[++cnt]=y;
	ne[cnt]=h[x];
	val[cnt]=z;
	val1[cnt]=a;
	h[x]=cnt;
	swap(x,y),z=0;
	v[++cnt]=y;
	ne[cnt]=h[x];
	val[cnt]=0;
	val1[cnt]=-a;
	h[x]=cnt;
}
void spfa(){
	memset(d,0x3f,sizeof d);
	queue<int>q;
	d[s]=0,f[s]=1;
	q.push(s);
	while(!q.empty()){
		long long u=q.front();
		q.pop();
		f[u]=0;
		for(long long i=h[u];i;i=ne[i]){
			if(val[i]==0)continue;
			long long vv=v[i],w=val1[i];
			if(d[vv]>d[u]+w){
				d[vv]=d[u]+w;
				if(!f[vv]){
					f[vv]=1;
					q.push(vv);
				}
			}
		}
	}
	for(int i=1;i<=n;i++)sn[i]=d[i];
}
struct st1{
	long long w,dd;
	st1(long long W,long long DD){
		w=W;
		dd=DD;
	}
	friend bool operator < (st1 x,st1 y){return x.dd>y.dd;}
};
bool bfs(){
	memset(d,0x3f,sizeof d);
	memset(f,0,sizeof f);
	for(int i=1;i<=n;i++)cur[i]=h[i];
	d[s]=0;
	priority_queue<st1>q;
	q.push(st1(s,0));
	while(q.size()){
		long long x=q.top().w;q.pop();
		if(f[x])continue;
		f[x]=1;
		for(long long i=h[x];i;i=ne[i]){
			if(val[i]==0)continue;
			long long vi=v[i],w=val1[i]+sn[x]-sn[vi];
			if(d[vi]>d[x]+w){
				d[vi]=d[x]+w;
				q.push(st1(vi,d[vi]));
			}
		}
	}
	for(long long i=1;i<=n;i++)d[i]=d[i]-sn[s]+sn[i];
	for(long long i=1;i<=n;i++)sn[i]=d[i];
	memset(f,0,sizeof f);
	return d[t]<0x3f3f3f3f3f3f;
}
long long dfs(long long x,long long y){
	if(x==t||y==0){
		return y;
	}
	f[x]=1;
	long long rp=0;
	for(long long i=cur[x];i;i=ne[i]){
		if(f[v[i]]==0&&d[v[i]]==d[x]+val1[i]&&val[i]){
			cur[x]=i;
			long long kk=dfs(v[i],min(y-rp,val[i]));
			rp+=kk;
			val[i]-=kk;
			val[i^1]+=kk;
			ans1+=kk*val1[i];
			if(y==rp)break;
		}
	}
	f[x]=0;
	return rp;
}
int main(){
	scanf("%lld%lld%lld%lld",&n,&m,&s,&t);
	for(long long i=1,u,v,w,x;i<=m;i++){
		scanf("%lld%lld%lld%lld",&u,&v,&w,&x);
		add(u,v,w,x);
	}
	spfa();
	while(bfs())ans+=dfs(s,0x3f3f3f3f3f3f3f3f);
	printf("%lld %lld",ans,ans1);
}
//注：在此代码中sn是势能函数，而h是链式前向星的head数组。