上文談到 網絡流-最大流 問題。
現在我們來學習 網絡流--費用流 這一塊,有紕漏的地方還請指出哦。
本文涉及的內容:
最大費用最大流
最小費用最大流
本文主要涉及的算法:
SPFA求費用流
dijkstra求費用流
zkw費用流
不明白最大流的讀者可以先去了解了解。
內容不會太難,畢竟作者能力有限,希望大家能有收獲。
費用流,也叫作最小費用最大流,是指在普通的網絡流圖中,每條邊的流量都有一個單價,求出一組可行解,使得在滿足它是最大流的情況下,總的費用最小。
首先開始最最小費用流:
最小費用最大流
顧名思義,就是利用最少的花費,達到流量最大的目的。
問題引入:
你是一個工廠的老板,達到可以制造任何多的貨物。眾所周知工廠不會見在市區,所以你需要運輸到市區,從工廠到銷售點有若干車次,每車次都有一個起點,終點(單向邊),還有一個容量。如今通常按勞分配,司機們取消了固定工資,取而代之的是運費,現在每輛車還有一個單位運費,指你運一單位貨物需要支付的錢,問當你滿足最大貨物運輸量的同時最少支出多少錢?

■ 每條邊多了一個值稱為費用。
■ 在最大化流量的同時最小化每條邊的費用與流量的乘積和。
■ 每次按照費用增廣?
■ 反向邊的費用設為原邊的相反數(想象司機運錯了給你再運回來並且會退錢 [ 好人吶!!])
■每次增廣的時候,流量$+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怎么辦?好說,滿足你。
■ $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; }
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\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 之內.。
zkw的慢
實踐中, 上面的這個算法非常奇怪. 在某一些圖上, 算法速度非常快, 另一些圖上卻比純 SPFA 增廣的算法慢. 不少同學經過實測總結的結果是稠密圖上比較快, 稀疏圖上比較慢, 但也不盡然. 這里我從理論上分析一下, 究竟這個算法用於哪些圖可以得到理想的效果.
先分析算法的增廣流程. 和 SPFA 直接算法相比, 由於同屬於沿最短路增廣的算法, 實際進行的增流操作並沒有太多的區別, 每次的增流路徑也大同小異. 因此不考慮多路增廣時, 增廣次數應該基本相同. 運行時間上主要的差異應當在於如何尋找增流路徑的部分.
那么 zkw 算法的優勢在於哪里呢? 與 SPFA 相比, KM 的重標號方式明顯在速度上占優, 每次只是一個對邊的掃描操作而已. 而 SPFA 需要維護較為復雜的標號和隊列操作, 同時為了修正標號, 需要不止一次地訪問某些節點, 速度會慢不少. 另外, 在 zkw 算法中, 增廣是多路進行的, 同時可能在一次重標號后進行多次增廣. 這個特點可以在許多路徑都費用相同的時候派上用場, 進一步減少了重標號的時間耗費.
下面想一想 zkw 算法的劣勢, 也就是 KM 重標號方式存在的問題. KM 重標號的主要問題就是, 不保證經過一次重標號之后能夠存在增廣路. 最差情況下, 一次只能在零權網絡中增加一條邊而已. 這時算法就會反復重標號, 反復嘗試增廣而次次不能增廣, 陷入弄巧成拙的境地.
接下來要說什么, 大家可能已經猜到了. 對於最終流量較大, 而費用取值范圍不大的圖, 或者是增廣路徑比較短的圖 (如二分圖), zkw 算法都會比較快. 原因是充分發揮優勢. 比如流多說明可以同一費用反復增廣, 費用窄說明不用改太多距離標號就會有新增廣路, 增廣路徑短可以顯著改善最壞情況, 因為即使每次就只增加一條邊也可以很快湊成最短路. 如果恰恰相反, 流量不大, 費用不小, 增廣路還較長, 就不適合 zkw 算法了.
本文部分材料來源於網絡,若有不當之處,會及時更改。
文章原創,轉載請標明出處,謝謝。
