上文談到 網絡流-最大流 問題。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 算法了.
本文部分材料來源於網絡,如有不當之處,會及時更改。
文章原創,轉載請標明出處,謝謝。