上文谈到 网络流-最大流 问题。html
如今咱们来学习 网络流--费用流 这一块,有纰漏的地方还请指出哦。ios
不明白最大流的读者能够先去了解了解。算法
内容不会太难,毕竟做者能力有限,但愿你们能有收获。编程
费用流,也叫做最小费用最大流,是指在普通的网络流图中,每条边的流量都有一个单价,求出一组可行解,使得在知足它是最大流的状况下,总的费用最小。网络
首先开始最最小费用流:ide
顾名思义,就是利用最少的花费,达到流量最大的目的。函数
你是一个工厂的老板,达到能够制造任何多的货物。众所周知工厂不会见在市区,因此你须要运输到市区,从工厂到销售点有若干车次,每车次都有一个起点,终点(单向边),还有一个容量。现在一般按劳分配,司机们取消了固定工资,取而代之的是运费,如今每辆车还有一个单位运费,指你运一单位货物须要支付的钱,问当你知足最大货物运输量的同时最少支出多少钱?学习
■ 每条边多了一个值称为费用。spa
■ 在最大化流量的同时最小化每条边的费用与流量的乘积和。code
■ 每次按照费用增广?
■ 反向边的费用设为原边的相反数(想象司机运错了给你再运回来而且会退钱 [ 好人呐!!])
■每次增广的时候,流量$+m$,那么费用增长$m \times 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)\n",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\n",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$不能处理负边权怎么办?
考虑给每一个点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\n",flow,res); return 0; }
最小费用流在 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\n",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 \notin V, uij > 0 \} \left \{ cij - Di + Dj } \right \}$ , 全部访问过的点距离标号增长
. 能够证实, 这样不会破坏性质 1, 并且至少有一条新的边进入了
的子图.
算法的步骤就是初始标号设为 , 不断增广, 若是不能增广, 修改标号继续增广, 直到完全不能增广: 源点的标号已经被加到了
. 注意: 在程序中全部的 cost 均表示的是 reduced cost, 即
. 另外, 这个算法不能直接用于有任何负权边的图. 更不能用于负权圈的状况. 有关这两种状况的处理, 参见 (2.) 和 (3.) 中的说明.
这样咱们获得了一个简单的算法, 只须要增广, 改标号, 各自只有 7 行, 不须要 BFS, 队列, SPFA, 编程复杂度很低. 因为实际的增广都是沿最短路进行的, 因此理论时间复杂度与使用 SPFA 等等方法的连续最短路算法一致, 但节省了 SPFA 或者 Dijkstra 的运算时间. 实测发现这种算法常数很小, 速度较快, employee 这道题全部数据加在一块儿耗时都在 2s 以内.。
实践中, 上面的这个算法很是奇怪. 在某一些图上, 算法速度很是快, 另外一些图上却比纯 SPFA 增广的算法慢. 很多同窗通过实测总结的结果是稠密图上比较快, 稀疏图上比较慢, 但也不尽然. 这里我从理论上分析一下, 究竟这个算法用于哪些图能够获得理想的效果.
先分析算法的增广流程. 和 SPFA 直接算法相比, 因为同属于沿最短路增广的算法, 实际进行的增流操做并无太多的区别, 每次的增流路径也大同小异. 所以不考虑多路增广时, 增广次数应该基本相同. 运行时间上主要的差别应当在于如何寻找增流路径的部分.
那么 zkw 算法的优点在于哪里呢? 与 SPFA 相比, KM 的重标号方式明显在速度上占优, 每次只是一个对边的扫描操做而已. 而 SPFA 须要维护较为复杂的标号和队列操做, 同时为了修正标号, 须要不止一次地访问某些节点, 速度会慢很多. 另外, 在 zkw 算法中, 增广是多路进行的, 同时可能在一次重标号后进行屡次增广. 这个特色能够在许多路径都费用相同的时候派上用场, 进一步减小了重标号的时间耗费.
下面想想 zkw 算法的劣势, 也就是 KM 重标号方式存在的问题. KM 重标号的主要问题就是, 不保证通过一次重标号以后可以存在增广路. 最差状况下, 一次只能在零权网络中增长一条边而已. 这时算法就会反复重标号, 反复尝试增广而次次不能增广, 陷入弄巧成拙的境地.
接下来要说什么, 你们可能已经猜到了. 对于最终流量较大, 而费用取值范围不大的图, 或者是增广路径比较短的图 (如二分图), zkw 算法都会比较快. 缘由是充分发挥优点. 好比流多说明能够同一费用反复增广, 费用窄说明不用改太多距离标号就会有新增广路, 增广路径短能够显著改善最坏状况, 由于即便每次就只增长一条边也能够很快凑成最短路. 若是偏偏相反, 流量不大, 费用不小, 增广路还较长, 就不适合 zkw 算法了.
本文部分材料来源于网络,如有不当之处,会及时更改。
文章原创,转载请标明出处,谢谢。