上文谈到 网络流-最大流 问题。
现在我们来学习 网络流--费用流 这一块，有纰漏的地方还请指出哦。
本文涉及的内容：
最大费用最大流
最小费用最大流
本文主要涉及的算法：
SPFA求费用流
dijkstra求费用流
zkw费用流
不明白最大流的读者可以先去了解了解。
内容不会太难，毕竟作者能力有限，希望大家能有收获。
费用流，也叫作最小费用最大流，是指在普通的网络流图中，每条边的流量都有一个单价，求出一组可行解，使得在满足它是最大流的情况下，总的费用最小。
首先开始最最小费用流：
最小费用最大流
顾名思义，就是利用最少的花费，达到流量最大的目的。
问题引入：
你是一个工厂的老板，达到可以制造任何多的货物。众所周知工厂不会见在市区，所以你需要运输到市区，从工厂到销售点有若干车次，每车次都有一个起点，终点(单向边)，还有一个容量。如今通常按劳分配，司机们取消了固定工资，取而代之的是运费，现在每辆车还有一个单位运费，指你运一单位货物需要支付的钱，问当你满足最大货物运输量的同时最少支出多少钱？
■ 每条边多了一个值称为费用。
■ 在最大化流量的同时最小化每条边的费用与流量的乘积和。
■ 每次按照费用增广？
■ 反向边的费用设为原边的相反数(想象司机运错了给你再运回来并且会退钱 [ 好人呐！！])
■每次增广的时候，流量$+m$，那么费用增加$m imes dis[ t ]$ (t为目标)
可行性
■ 会不会出现负环？
只要初始时没有负环，之后就不会有负环，因为真想和反向花费相反。
■ 注意到初始时所有反向边的残量为0，可以只考虑原图中的边。
■ 如果增广路中出现了负环，那么在上一次选择中一定有一条更短的路径。
■ 如果开始就有负环呢？
那么它说明你图建错了(存在某种玄学的消环的方法，但是好麻烦QAQ，而且时间复杂度得不到保证。)
正确性
■ 为什么是对的？
当前最小费用流 <=> 残量网络无负环
如果有负环我们可以从这个负环上走一圈来得到更小的解。
这东西反过来也是成立的，即如果有更小解，一定存在一个负环来让我们走一圈。
时间复杂度
这我咋知道，网络流的时间复杂度只有$O(TLE)$和$O(not TLE)$
一般的网络流根本跑不到上界，如果想要了解跑到上界的算法，可以去了解“前置-重贴标签算法”。这里从网上找了一份代码：
#include <stdio.h> #include <stdlib.h> #include <limits.h> //图节点 typedef struct AdjacentNodeType { int index; struct AdjacentNodeType *nextNeighbor; }AdjacentNode,*pAdjacentNode; typedef struct VertexNode { int index; char name; int h; int e; pAdjacentNode current; pAdjacentNode head; struct VertexNode *next; struct VertexNode *pre; }Vertex,*pVertex; //图 typedef struct { int vn; int **E; //容量C作为边矩阵的值 pVertex *V; int **f; //流值 int **rE; // 残留边 pVertex L; //前置重贴标签用到的链表，在initGraph()中初始化 }Graph,*pGraph; void generateAdjacentList(pGraph g,int s,int t) { for(int i=0;i<g->vn;i++) { for(int j=0;j<g->vn;j++) { if(g->E[i][j]>0 || g->E[j][i]>0) { pAdjacentNode adj=(pAdjacentNode)malloc(sizeof(AdjacentNode)); adj->index=j; adj->nextNeighbor=g->V[i]->head; g->V[i]->head=adj; } } } } //根据算法导论 图26-6初始化图 pGraph initGraph() { pGraph g=(pGraph)malloc(sizeof(Graph)); g->vn=6; int s=0,t=g->vn-1; pVertex vs=(pVertex)malloc(sizeof(Vertex)); vs->name='s'; vs->index=0; pVertex v1=(pVertex)malloc(sizeof(Vertex)); v1->name='1'; v1->index=1; pVertex v2=(pVertex)malloc(sizeof(Vertex)); v2->name='2'; v2->index=2; pVertex v3=(pVertex)malloc(sizeof(Vertex)); v3->name='3'; v3->index=3; pVertex v4=(pVertex)malloc(sizeof(Vertex)); v4->name='4'; v4->index=4; pVertex vt=(pVertex)malloc(sizeof(Vertex)); vt->name='t'; vt->index=5; g->V=(pVertex*)malloc(g->vn*sizeof(pVertex)); g->V[0]=vs; g->V[1]=v1; g->V[2]=v2; g->V[3]=v3; g->V[4]=v4; g->V[5]=vt; for(int i=0;i<g->vn;i++) { g->V[i]->current=g->V[i]->head=NULL; g->V[i]->pre=g->V[i]->next=NULL; } g->L=NULL; for(int i=g->vn-1;i>=0;i--) { if(i==s || i==t) continue; g->V[i]->next=g->L; if(g->L) g->L->pre=g->V[i]; g->L=g->V[i]; } g->E = (int**)malloc(g->vn*sizeof(int*)); g->f= (int**)malloc(g->vn*sizeof(int*)); g->rE= (int**)malloc(g->vn*sizeof(int*)); for(int i=0;i<g->vn;i++) { g->E[i]=(int*)malloc(g->vn*sizeof(int)); g->f[i]=(int*)malloc(g->vn*sizeof(int)); g->rE[i]=(int*)malloc(g->vn*sizeof(int)); } for(int i=0;i<g->vn;i++) { for(int j=0;j<g->vn;j++) { g->E[i][j]=0; //g->f[i][j]=0; } } g->E[0][1]=16; g->E[0][2]=13; g->E[1][3]=12; g->E[2][1]=4; g->E[2][4]=14; g->E[3][2]=9; g->E[3][5]=20; g->E[4][3]=7; g->E[4][5]=4; generateAdjacentList(g,s,t); return g; } void initResidualNetwork(pGraph g) { for(int i=0;i<g->vn;i++) { for(int j=0;j<g->vn;j++) { g->rE[i][j]=0; } } for(int i=0;i<g->vn;i++) { for(int j=0;j<g->vn;j++) { if(g->E[i][j]>0) { int x=g->E[i][j]-g->f[i][j]; if(x>0) g->rE[i][j]=x; if(g->f[i][j]>0) g->rE[j][i]=g->f[i][j]; } } } } void initializePreflow(pGraph g,int s) { for(int i=0;i<g->vn;i++) { g->V[i]->e=g->V[i]->h=0; } for(int i=0;i<g->vn;i++) { for(int j=0;j<g->vn;j++) { g->f[i][j]=0; } } g->V[s]->h=g->vn; for(int i=0;i<g->vn;i++) { if(i != s) { g->f[s][i]=g->E[s][i]; g->V[i]->e=g->E[s][i]; g->V[s]->e -= g->E[s][i]; } } } void push(pGraph g,int u,int v) { if(g->V[u]->e<=0) return; if(g->V[u]->h != g->V[v]->h+1) return; if(g->rE[u][v]<=0) return; int d=g->V[u]->e < g->rE[u][v] ? g->V[u]->e : g->rE[u][v]; //更新流 if(g->E[u][v]>0) { g->f[u][v]+=d; } else { g->f[v][u]-=d; } //更新超额流 g->V[u]->e-=d; g->V[v]->e+=d; //更新残存图 if(g->E[u][v]>0) { int x=g->E[u][v]-g->f[u][v]; g->rE[u][v]=x; if(g->f[u][v]>0) g->rE[v][u]=g->f[u][v]; } } //进入函数时，默认保证g->V[u]->e>0 //返回能从u进行push的邻接顶点位置 //返回-1代表残留图中该顶点无邻接点 int relabel(pGraph g,int u) { int minh=INT_MAX,minPos; for(int i=0;i<g->vn;i++) { if(i==u) continue; if(g->rE[u][i]>0) { if(g->V[i]->h<g->V[u]->h) return i; else { if(g->V[i]->h<minh) { minh=g->V[i]->h; minPos=i; } } } } if(minh<INT_MAX) { g->V[u]->h=minh+1; return minPos; } else//u没有邻接点时走到这里 return -1; } void discharge(pGraph g,int u) { while(g->V[u]->e>0) { pAdjacentNode pv=g->V[u]->current; if(pv==NULL) { relabel(g,u); g->V[u]->current=g->V[u]->head; } else if(g->rE[u][pv->index]>0 && g->V[u]->h == g->V[pv->index]->h+1) { push(g,u,pv->index); } else g->V[u]->current=pv->nextNeighbor; } } void relabelToFront(pGraph g,int s,int t) { initializePreflow(g,s); initResidualNetwork(g); for(int i=0;i<g->vn;i++) { if(i==s || i==t) continue; g->V[i]->current=g->V[i]->head; } pVertex pu=g->L; while(pu!=NULL) { int oldHeight=pu->h; discharge(g,pu->index); if(pu->h>oldHeight) { //链表中需要移动的节点是头节点则不移动 if(pu!=g->L) { pu->pre->next=pu->next; if(pu->next) pu->next->pre=pu->pre; pu->next=g->L; g->L->pre=pu; pu->pre=NULL; g->L=pu; } } pu=pu->next; } } void printFlow(pGraph g) { for(int i=0;i<g->vn;i++) { for(int j=0;j<g->vn;j++) { if(g->E[i][j]>0) printf("%c->%c (%d/%d) ",g->V[i]->name,g->V[j]->name,g->f[i][j],g->E[i][j]); } } } int calcuMaxFlow(pGraph g,int s) { int maxFlow=0; for(int i=0;i<g->vn;i++) { if(i==s) continue; if(g->f[s][i]>0) maxFlow+=g->f[s][i]; } return maxFlow; } void main() { pGraph g=initGraph(); int s=0,t=g->vn-1; relabelToFront(g,s,t); int maxFlow=calcuMaxFlow(g,s); printf("Max Flow is:%d ",maxFlow); printFlow(g); getchar(); }
实现
我们考虑一下我们是怎么做最大流的，我们是将增广路按照距离来$bfs$分层，那么这个我们也可以模仿此，但是每次我们怎么走呢？
spfa中，我们按照费用的最小来走，这样子的话，就很明显是为了走最小花费了。
每次不断的去找最短路，从后往前依次更新用到的边。
最大流中已经提到了解法，在这里就不过多解释了。
$bfs$变成了$spfa$返回值还是一样的。
伪代码：
int MCMF(int S, int T) { int ans = 0; while (SPFA(S, T)) { int minv = 0x7fffffff; for (int x = T; x != S; x = from[fa[x]]) minv = min(minv, ret[fa[x]]); ans += minv * dis[T]; for (int x = T; x != S; x = from[fa[x]]) { ret[fa[x]] -= minv; ret[rev[fa[x]]] += minv; } } return ans; }
真●代码：
bool bfs(int s,int t) { memset(flow,0x7f,sizeof(flow)); memset(dis,0x7f,sizeof(dis)); memset(vis,0,sizeof(vis)); Q.push(s);vis[s]=1;dis[s]=0,pre[t]=-1; while(!Q.empty()) { int temp=Q.front(); Q.pop(); vis[temp]=0; for(int i=head[temp];i!=-1;i=edge[i].nxt) { int v=edge[i].to; if(edge[i].flow>0&&dis[v]>dis[temp]+edge[i].dis) { dis[v]=dis[temp]+edge[i].dis; pre[v]=temp; last[v]=i; flow[v]=min(flow[temp],edge[i].flow); if(!vis[v]) { vis[v]=1; Q.push(v); } } } } return pre[t]!=-1; } void MCMF() { while(bfs(s,t)) { int now=t; maxflow+=flow[t]; mincost+=flow[t]*dis[t]; while(now!=s) { edge[last[now]].flow-=flow[t]; edge[last[now]^1].flow+=flow[t]; now=pre[now]; } } }
想写dijkstra怎么办？好说，满足你。
■ $dijkstra$不能处理负边权怎么办？
考虑给每个点x,加一个“势”hx。一条u -> v的费用为 c 的边，变成一条u -> v费用是&c-hx+hu&的边。
■$a -> b -> c ->... -> z$的路径增加了$(ha-hb)+(hb-hc)+....+(hy-hz) = ha-hz$.
这告诉我们，这样处理的最短路和原来的最短路是相同的只是增加了$ha-hz$
■ 如果我们增广时能给每个点找到一个势，使得所有边处理之后费用非负，那么就可以用$dijkstra$来求最短路了。
$$c − hv + hu ≥ 0 ⇔ hv ≤ c + hu$$
上边这个式子看上去很像一个单源最短路。但我们不能直接使用单源最短路的值，因为我们正要求它。
能不能使用上一次最短路的值？
答案是可以，也就是说我们每次增广的算法流程就是：
1.领本次的势hx为上一次的disx
2.利用带势函数的$dijkstra$求出最短路。
3,.增广这条最短路。
正确性
这为什么是对的呢？
考虑一条边 u → v ，费用为 c 。
如果它上一次增广时残量不为 0 ，那么根据最短路的性质有
$dis[u] + c >= dis[v]$
不然说明最短求错了。
如果它上次增⼴时残量为 0 而现在不为 0 ，那说明它的反向边被增广了。而增广的路径是最短路径，反向边是 v → u ，费用为 −c 。所以$dis[u] = dis[v] - c$，也就是说$c + dis[u] − dis [v] = 0$ ，也是非负的。
于是乎我们就可以用$dijkstra$增广最短路了，而且很难卡。
$dijkstra$时每条边的长度看作其费用$+dis[u]-dis[v]$. $dijkstra$结束前将$dis[x]+= hx$
int MCMF(int S, int T) { int ans = 0; while (SPFA(S, T)) { int minv = 0x7fffffff; for (int x = T; x != S; x = from[fa[x]]) minv = min(minv, ret[fa[x]]); ans += minv * dis[T]; for (int x = T; x != S; x = from[fa[x]]) { ret[fa[x]] -= minv; ret[rev[fa[x]]] += minv; } for (int x = 1; x <= n; x++) h[x] = dis[x]; } return ans; }
#pragma GCC optimize(2) #include <iostream> #include <algorithm> #include <queue> #include <math.h> #include <stdio.h> #include <string.h> #include <algorithm> #define MAXN_ 5050 #define INF 0x3f3f3f3f #define P pair<int,int> using namespace std; struct edge { int to,cap,cost,rev; }; int n,m,flow,s,t,cap,res,cost,from,to,h[MAXN_]; std::vector<edge> G[MAXN_]; int dist[MAXN_],prevv[MAXN_],preve[MAXN_]; // 前驱节点和对应边 inline void add() { G[from].push_back((edge) { to,cap,cost,(int)G[to].size() }); G[to].push_back((edge) { from,0,-cost,(int)G[from].size()-1 }); } // 在vector 之中找到边的位置所在! inline int read() { int x=0; char c=getchar(); bool flag=0; while(c<'0'||c>'9') { if(c=='-')flag=1; c=getchar(); } while(c>='0'&&c<='9') { x=(x<<3)+(x<<1)+c-'0'; c=getchar(); } return flag?-x:x; } inline void min_cost_flow(int s,int t,int f) { fill(h+1,h+1+n,0); while(f > 0) { priority_queue<P,vector<P>, greater<P> > D; memset(dist,INF,sizeof dist); dist[s] = 0; D.push(P(0,s)); while(!D.empty()) { P now = D.top(); D.pop(); if(dist[now.second] < now.first) continue; int v = now.second; for(int i=0; i<(int)G[v].size(); ++i) { edge &e = G[v][i]; if(e.cap > 0 && dist[e.to] > dist[v] + e.cost + h[v] - h[e.to]) { dist[e.to] = dist[v] + e.cost + h[v] - h[e.to]; prevv[e.to] = v; preve[e.to] = i; D.push(P(dist[e.to],e.to)); } } } // 无法增广 ， 就是找到了答案了！ if(dist[t] == INF) break; for(int i=1; i<=n; ++i) h[i] += dist[i]; int d = f; for(int v = t; v != s; v = prevv[v]) d = min(d,G[prevv[v]][preve[v]].cap); f -= d; flow += d; res += d * h[t]; for(int v=t; v!=s; v=prevv[v]) { edge &e = G[prevv[v]][preve[v]]; e.cap -= d; G[v][e.rev].cap += d; } } } int main(int argc,char const* argv[]) { n = read(); m = read(); s = read(); t = read(); for(int i=1; i<=m; ++i) { from = read(); to = read(); cap = read(); cost = read(); add(); } min_cost_flow(s,t,INF); printf("%d %d ",flow,res); return 0; }
zkw算法[转]
最小费用流在 OI 竞赛中应当算是比较偏门的内容, 但是 NOI2008 中 employee 的突然出现确实让许多人包括 zkw 自己措手不及. 可怜的 zkw 当时想出了最小费用流模型, 可是他从来没有实现过, 所以不敢写, 此题 0 分. zkw 现在对费用流的心得是: 虽然理论上难, 但是写一个能 AC 题的费用流还算简单. 先贴一个我写的 costflow 程序: 只有不到 70 行, 费用流比最大流还好写～。
#include <cstdio> #include <cstring> using namespace std; const int maxint=~0U>>1; int n,m,pi1,cost=0; bool v[550]; struct etype { int t,c,u; etype *next,*pair; etype() {} etype(int t_,int c_,int u_,etype* next_): t(t_),c(c_),u(u_),next(next_) {} void* operator new(unsigned,void* p) { return p; } } *e[550]; int aug(int no,int m) { if(no==n)return cost+=pi1*m,m; v[no]=true; int l=m; for(etype *i=e[no]; i; i=i->next) if(i->u && !i->c && !v[i->t]) { int d=aug(i->t,l<i->u?l:i->u); i->u-=d,i->pair->u+=d,l-=d; if(!l)return m; } return m-l; } bool modlabel() { int d=maxint; for(int i=1; i<=n; ++i)if(v[i]) for(etype *j=e[i]; j; j=j->next) if(j->u && !v[j->t] && j->c<d)d=j->c; if(d==maxint)return false; for(int i=1; i<=n; ++i)if(v[i]) for(etype *j=e[i]; j; j=j->next) j->c-=d,j->pair->c+=d; pi1 += d; return true; } int main() { freopen("costflow.in","r",stdin); freopen("costflow.out","w",stdout); scanf("%d %d",&n,&m); etype *Pe=new etype[m+m]; while(m--) { int s,t,c,u; scanf("%d%d%d%d",&s,&t,&u,&c); e[s]=new(Pe++)etype(t, c,u,e[s]); e[t]=new(Pe++)etype(s,-c,0,e[t]); e[s]->pair=e[t]; e[t]->pair=e[s]; } do do memset(v,0,sizeof(v)); while(aug(1,maxint)); while(modlabel()); printf("%d ",cost); return 0; }
这里使用的是连续最短路算法. 最短路算法? 为什么程序里没有 SPFA? Dijkstra? 且慢, 先让我们回顾一下图论中最短路算法中的距离标号. 定义 
为点 
的距离标号, 任何一个最短路算法保证, 算法结束时对任意指向顶点 、从顶点 
出发的边满足 
(条件1), 且对于每个 存在一个 使得等号成立 (条件2). 换句话说, 任何一个满足以上两个条件的算法都可以叫做最短路, 而不仅仅是 SPFA、Dijkstra, 算法结束后, 恰在最短路上的边满足 
.
在最小费用流的计算中, 我们每次沿 的路径增广后都不会破坏条件 1, 但是可能破坏了条件 2. 不满足条件 2 的后果是什么呢? 使我们找不到每条边都满足 
新的增广路. 只好每次增广后使用 Dijkstra, SPFA 等等算法重新计算新的满足条件 2 的距离标号. 这无疑是一种浪费. KM 算法中我们可以修改不断修改可行顶标, 不断扩大可行子图, 这里也同样, 我们可以在始终满足条件 1 的距离标号上不断修改, 直到可以继续增广 (满足条件 2).
回顾一下 KM 算法修改顶标的方法. 根据最后一次寻找交错路不成功的 DFS, 找到 
, 左边的点增加 
, 右边的点减少 . 这里也一样, 根据最后一次寻找增广路不成功的 DFS, 找到 $ d = min {i in V, j
otin V, uij > 0 } left { cij - Di + Dj }
ight }$ , 所有访问过的点距离标号增加 . 可以证明, 这样不会破坏性质 1, 而且至少有一条新的边进入了 的子图.
算法的步骤就是初始标号设为 
, 不断增广, 如果不能增广, 修改标号继续增广, 直到彻底不能增广: 源点的标号已经被加到了 
. 注意: 在程序中所有的 cost 均表示的是 reduced cost, 即 
. 另外, 这个算法不能直接用于有任何负权边的图. 更不能用于负权圈的情况. 有关这两种情况的处理, 参见 (2.) 和 (3.) 中的说明.
这样我们得到了一个简单的算法, 只需要增广, 改标号, 各自只有 7 行, 不需要 BFS, 队列, SPFA, 编程复杂度很低. 由于实际的增广都是沿最短路进行的, 所以理论时间复杂度与使用 SPFA 等等方法的连续最短路算法一致, 但节省了 SPFA 或者 Dijkstra 的运算时间. 实测发现这种算法常数很小, 速度较快, employee 这道题所有数据加在一起耗时都在 2s 之内.。
zkw的慢
实践中, 上面的这个算法非常奇怪. 在某一些图上, 算法速度非常快, 另一些图上却比纯 SPFA 增广的算法慢. 不少同学经过实测总结的结果是稠密图上比较快, 稀疏图上比较慢, 但也不尽然. 这里我从理论上分析一下, 究竟这个算法用于哪些图可以得到理想的效果.
先分析算法的增广流程. 和 SPFA 直接算法相比, 由于同属于沿最短路增广的算法, 实际进行的增流操作并没有太多的区别, 每次的增流路径也大同小异. 因此不考虑多路增广时, 增广次数应该基本相同. 运行时间上主要的差异应当在于如何寻找增流路径的部分.
那么 zkw 算法的优势在于哪里呢? 与 SPFA 相比, KM 的重标号方式明显在速度上占优, 每次只是一个对边的扫描操作而已. 而 SPFA 需要维护较为复杂的标号和队列操作, 同时为了修正标号, 需要不止一次地访问某些节点, 速度会慢不少. 另外, 在 zkw 算法中, 增广是多路进行的, 同时可能在一次重标号后进行多次增广. 这个特点可以在许多路径都费用相同的时候派上用场, 进一步减少了重标号的时间耗费.
下面想一想 zkw 算法的劣势, 也就是 KM 重标号方式存在的问题. KM 重标号的主要问题就是, 不保证经过一次重标号之后能够存在增广路. 最差情况下, 一次只能在零权网络中增加一条边而已. 这时算法就会反复重标号, 反复尝试增广而次次不能增广, 陷入弄巧成拙的境地.
接下来要说什么, 大家可能已经猜到了. 对于最终流量较大, 而费用取值范围不大的图, 或者是增广路径比较短的图 (如二分图), zkw 算法都会比较快. 原因是充分发挥优势. 比如流多说明可以同一费用反复增广, 费用窄说明不用改太多距离标号就会有新增广路, 增广路径短可以显著改善最坏情况, 因为即使每次就只增加一条边也可以很快凑成最短路. 如果恰恰相反, 流量不大, 费用不小, 增广路还较长, 就不适合 zkw 算法了.
本文部分材料来源于网络，若有不当之处，会及时更改。
文章原创，转载请标明出处，谢谢。