下面的最小费用最大流算法采用的是“基于Floyd最短路算法的Ford和Fulkerson迭加算法”,其基本思路为:把各条弧上单位流量的费用看成某种长度,用Floyd求最短路的方法确定一条自V1至Vn的最短路;再将这条最短路作为可扩充路,用求解最大流问题的方法将其上的流量增至最大可能值;而这条最短路上的流量增加后,其上各条弧的单位流量的费用要重新确定,如此多次迭代,最终得到最小费用最大流。. k6 X0 r( j6 f7 B c
! L1 a4 Y% C4 r8 r8 Jfunction [f,MinCost,MaxFlow]=MinimumCostFlow(a,c,V,s,t). Q8 L0 ~. U. X+ A/ B$ }
%% 基于Floyd最短路算法的Ford和Fulkerson迭加算法3 n2 f5 H, c7 l5 Q1 J
%% 输入参数列表 $ v Y& a9 q( R2 j' F% a 单位流量的费用矩阵 ! ?. G) J1 d) ?% c 链路容量矩阵 , T9 N6 h4 B) M7 t" d) S% V 最大流的预设值,可为无穷大& {1 p2 v; c+ ]/ T6 m$ S1 t
% s 源节点 " z. z, P' ]. A, {7 y% t 目的节点8 W; s1 k7 P1 G
%% 输出参数列表3 k7 S7 i$ d% w/ h
% f 链路流量矩阵; S2 V, V6 i" u5 C* \
% MinCost 最小费用 8 ~* B0 U. c4 G% MaxFlow 最大流量% j' T. k$ l7 r0 @" |
%% 第一步:初始化( T. t. f0 w- p
N=size(a,1);%节点数目) C$ u) T' I0 P/ t, X* y( `
f=zeros(N,N);%流量矩阵,初始时为零流 T& n# g$ D2 g7 q% ]' M
MaxFlow=sum(f(s,);%最大流量,初始时也为零 o& F% Y" H8 N+ a$ P
flag=zeros(N,N);%真实的前向边应该被记住. K' i6 ?+ t% z
for i=1:N, ^: @ n. X& _( j0 Q2 k9 D! k% b
for j=1:N - F& V& M% S$ G C; `3 q+ j. F6 Uif i~=j&&c(i,j)~=0% C/ T7 S# P# i. J! ], ~
flag(i,j)=1;%前向边标记 % f7 [5 ^/ q0 q8 X' q, _8 W- `, R% Fflag(j,i)=-1;%反向边标记, ~3 |( }* k1 e# r2 L8 F3 n' @
end $ p9 |6 {3 J" \* {; eif a(i,j)==inf 5 V8 Z6 P) b& Ha(i,j)=BV;- Y8 R/ b% D; {7 T U; b, t
w(i,j)=BV;%为提高程序的稳健性,以一个有限大数取代无穷大4 [' S/ \5 i% ?
end c' a4 J( d" w! c$ o F9 O* ` Zend 0 G# U; Z, Y. M/ oend * I3 |, y* t' }) R$ L. ]if L(end) RE=1;%如果路径长度小于大数,说明路径存在2 f; w# M$ x) {! m* ^# \: G
else5 ^: L# o8 F$ g" Y8 j
RE=0; 3 i9 c- _. Q: e0 Z* Z3 u5 ~$ Kend- @: R; G7 i: ?- ^$ V0 W( x; [9 @7 H: \+ U
%% 第二步:迭代过程 . i5 f3 R1 t D7 Zwhile RE==1&&MaxFlow<=V%停止条件为达到最大流的预设值或者没有从s到t的最短路: B* u/ i; v% E1 y# P
%以下为更新网络结构 9 ]0 D. h9 B* q! t5 R8 P7 @. ^MinCost1=sum(sum(f.*a));0 E' U: I% W, }" U5 G" p
MaxFlow1=sum(f(s,); 4 I$ T: r' q: K, W! v m0 kf1=f; ' \; r K( `. x8 @& ~TS=length(R)-1;%路径经过的跳数 d: b8 ]! p$ R. c: y* iLY=zeros(1,TS);%流量裕度& f4 C i. e# i' h" d
for i=1:TS! z" @. @! k, E9 u5 \3 @
LY(i)=c(R(i),R(i+1));6 `: x4 x( e b* A9 A
end - M2 S5 M- Q' t, ]$ `( OmaxLY=min(LY);%流量裕度的最小值,也即最大能够增加的流量 4 V$ T0 l. _1 S2 ~! Q# S/ L: M4 ~6 `for i=1:TS! V3 C' _8 S7 C7 I& N( u E! S
u=R(i); `" l$ d$ I) h( B7 Qv=R(i+1); ' A( t' w, k# l; eif flag(u,v)==1&&maxLY f(u,v)=f(u,v)+maxLY;%记录流量值 3 `) @+ Q5 v0 J# J9 N# _) T9 ww(u,v)=a(u,v);%更新权重值; W0 D4 a$ F4 Y/ W" h
c(v,u)=c(v,u)+maxLY;%反向链路的流量裕度更新 2 N# l5 _: A; I3 e9 z6 y5 G, r6 ?elseif flag(u,v)==1&&maxLY==c(u,v)%当这条边为前向边且是饱和边时4 n9 z, }- x8 _3 i; ]! s
w(u,v)=BV;%更新权重值9 D3 |, J7 b- k* J9 H) ^
c(u,v)=c(u,v)-maxLY;%更新流量裕度值 |4 q7 q# `& I) v, yw(v,u)=-a(u,v);%反向链路权重更新 0 P! j+ K5 I3 e4 Xelseif flag(u,v)==-1&&maxLY w(v,u)=a(v,u);) w; p+ e9 p" t7 ~
c(v,u)=c(v,u)+maxLY; 6 i1 a m- X5 j( u+ Jw(u,v)=-a(v,u); , v9 i1 d9 F7 H) y8 M2 ]% Ielseif flag(u,v)==-1&&maxLY==c(u,v)%当这条边为反向边且是饱和边时 5 }+ ~5 P" {% z2 j3 Z- u* E. Vw(v,u)=a(v,u); 2 J) {9 S& l. s6 Dc(u,v)=c(u,v)-maxLY;# S" R8 Q' f9 `7 B; h1 g
w(u,v)=BV; 9 v# W; R f7 ^- K& C2 a! [else- G5 a' @9 g( d) h. w4 j! R+ ?! h! T
end0 U; p% ]5 b: X# b
end 5 v( I5 B8 a8 _ WMaxFlow2=sum(f(s,);. C; A# o( X( o* b& R
MinCost2=sum(sum(f.*a)); 7 a& ~& w- Z8 P# k! z pif MaxFlow2<=V 1 y# n$ j4 s7 NMaxFlow=MaxFlow2; ; l( o. T2 `3 Y2 j2 ]MinCost=MinCost2;: W! t9 t& u$ R, u4 K3 V
[L,R]=FLOYD(w,s,t);7 j0 f5 O1 I$ A0 ?
else w$ `( p \6 C0 G$ A
f=f1+prop*(f-f1); ' N3 z- C- F6 @0 e+ PMaxFlow=V; " |4 O5 l; s+ ~4 U$ cMinCost=MinCost1+prop*(MinCost2-MinCost1); - Y3 G/ H4 q- B" C2 Z- Treturn; V) A9 T) B" L8 L4 v4 H4 Y4 N+ c
end + P/ I7 ^' H; U8 f4 h9 ~if L(end) RE=1;%如果路径长度小于大数,说明路径存在$ W3 m" o6 v1 A% S4 J
else 4 t* ]8 M; ^+ hRE=0; ! i# C$ S* {+ J* V5 K3 ^. K7 b+ wend ( `, W+ f8 }7 Mend 0 @) ^" X1 x9 ?function [L,R]=FLOYD(w,s,t) # u+ Q4 l$ j+ y- O* V! m+ U5 N, zn=size(w,1);5 i8 ]6 z% c. A ~7 l! Z
D=w; ( ^" u* g |. q" m8 P1 |( _& q Wpath=zeros(n,n);5 w) f3 j! G" x+ f
%以下是标准floyd算法3 u: K$ Y) l6 t$ G3 A
for i=1:n - @1 `6 s* H- g; k2 Qfor j=1:n : k% ]# L& }4 e" z+ {' ~2 `6 Pif D(i,j)~=inf k. ?/ r; \ v5 W7 _3 t0 k
path(i,j)=j; 4 K; D3 H. c# j; Gend ) }: n# |8 G& O1 R( J- {# ?" xend3 }) J* A5 {. ]! [, I; V+ M1 L
end3 L& W/ ^- x$ r" L
for k=1:n# R, _) H, o8 ~3 V8 M1 H
for i=1:n # ~$ F- N8 p) R8 s5 ?) Pfor j=1:n : y) W( h$ R% I. k( Cif D(i,k)+D(k,j) D(i,j)=D(i,k)+D(k,j);) E! f. a* N3 T/ @$ ]
path(i,j)=path(i,k);& G+ K* a+ }2 h: `, j/ U# G
end7 h# @( [8 A+ u3 R8 t1 i% u0 I
end6 j. g% L3 Y4 U2 C$ p+ H
end $ t* s, r; L: Q2 cend# G% ]/ S; P: n) m) S. ]
L=zeros(0,0);3 ]7 E9 t- A$ K# G6 ~& Y
R=s;6 w& s. f2 P$ F3 Y( W% `! h
while 1 * T4 {8 F2 K: e" S) O8 hif s==t: g# s, G/ x* l: }
L=fliplr(L);- Q' w7 v/ ]. Z) ?
L=[0,L]; & d; b k) e |; kreturn q/ H6 @4 L. R G, y8 i( M5 p7 Lend 0 m, h6 I9 L$ r& u) z4 W) iL=[L,D(s,t)];$ Q: Y& I O5 F* h8 a' a
R=[R,path(s,t)]; ) ^2 V" M/ g1 z- v7 D6 [8 I4 {s=path(s,t);5 A6 M# z9 J t0 r/ x
end