下面的最小费用最大流算法采用的是“基于Floyd最短路算法的Ford和Fulkerson迭加算法”,其基本思路为:把各条弧上单位流量的费用看成某种长度,用Floyd求最短路的方法确定一条自V1至Vn的最短路;再将这条最短路作为可扩充路,用求解最大流问题的方法将其上的流量增至最大可能值;而这条最短路上的流量增加后,其上各条弧的单位流量的费用要重新确定,如此多次迭代,最终得到最小费用最大流。; n/ ^; J. ~- _1 i
9 {7 Y% M; C5 @( B8 Nfunction [f,MinCost,MaxFlow]=MinimumCostFlow(a,c,V,s,t) 4 b! S& W9 G& d% D%% 基于Floyd最短路算法的Ford和Fulkerson迭加算法 7 `$ a, x0 [0 J7 V% I6 _" w6 _%% 输入参数列表$ v5 V! [; }2 R+ l9 }
% a 单位流量的费用矩阵4 @" t: J3 r7 Z: v- K9 {( Y
% c 链路容量矩阵, q: u7 }' t3 _1 [( z5 c
% V 最大流的预设值,可为无穷大4 w R! s. M6 H
% s 源节点 ! g: O$ _. u; }% t 目的节点 - ^2 d8 `" X" \2 {0 h3 X) e%% 输出参数列表& _' Q3 Y* O1 T
% f 链路流量矩阵 6 _5 J* B2 o2 T0 o) b% MinCost 最小费用. s- ?7 Y1 }& O# r) p' l
% MaxFlow 最大流量 3 \) Y/ C- [$ D" F%% 第一步:初始化9 L. y) y8 M0 C. T p9 {- c
N=size(a,1);%节点数目5 O/ \/ C% e) T. W
f=zeros(N,N);%流量矩阵,初始时为零流. x$ e8 _1 W3 N4 N; [
MaxFlow=sum(f(s,);%最大流量,初始时也为零. x6 Y5 {2 N1 b4 j( ~( x
flag=zeros(N,N);%真实的前向边应该被记住 7 m, @/ Y& E' R# Sfor i=1:N. g7 N/ `: A4 Q2 Y* O# m
for j=1:N7 Q; p9 W9 o2 D! \! q2 J
if i~=j&&c(i,j)~=0 6 ]( ^& Z% j3 ~: ]1 k2 @/ Iflag(i,j)=1;%前向边标记& S l& Y* N! E/ U/ `0 f
flag(j,i)=-1;%反向边标记8 N) [% v+ X7 e3 y
end! }' _' A1 ]5 J6 ]# Q2 Y
if a(i,j)==inf 6 B8 z1 b8 y7 D' Ba(i,j)=BV; * B7 w) o$ M) \3 \* M+ A# p, ]w(i,j)=BV;%为提高程序的稳健性,以一个有限大数取代无穷大 3 K/ ?' ^4 W" V" v7 r5 m) G3 eend b7 h. M0 I5 H# g! pend ; ?! |% o: I- t2 Send( _ Q2 v/ g0 l) r
if L(end) RE=1;%如果路径长度小于大数,说明路径存在 2 Y0 I8 c* u7 W- kelse# ]7 H6 f6 G9 C. R/ l
RE=0; ; k7 U- o% C2 I1 Q" N; Xend L5 @/ X9 s# z7 F+ P1 V
%% 第二步:迭代过程 " ` a/ R5 h7 ]8 ~6 w- E& Swhile RE==1&&MaxFlow<=V%停止条件为达到最大流的预设值或者没有从s到t的最短路 4 N3 ]3 L" t! x' B1 D: _5 V* N4 i%以下为更新网络结构2 A* I @+ |! l% h# G2 m& x
MinCost1=sum(sum(f.*a));! W% b0 ~+ Z; m0 K2 X# ^
MaxFlow1=sum(f(s,);) s; i' ~- ~. O' m2 [1 x7 Q
f1=f; 9 p! O0 w w, i/ F9 w3 |$ GTS=length(R)-1;%路径经过的跳数* L! _: ?! ?) C/ }+ ^9 g0 p! ~
LY=zeros(1,TS);%流量裕度 4 I) g+ D7 g* V* ?for i=1:TS $ G8 W( I) H8 ?, lLY(i)=c(R(i),R(i+1));# m% k5 Z7 M6 S; m: H) E
end 6 k4 }5 f' d2 `1 K7 |: ImaxLY=min(LY);%流量裕度的最小值,也即最大能够增加的流量* `% Y! i# q/ E' b6 o7 a& D+ K
for i=1:TS# c+ c3 Y8 Y P
u=R(i);9 Z) |& \/ }! _7 z
v=R(i+1);6 z- K l6 M/ }
if flag(u,v)==1&&maxLY f(u,v)=f(u,v)+maxLY;%记录流量值$ w; R5 _$ |* k
w(u,v)=a(u,v);%更新权重值 }2 W/ _/ e9 j- `1 _& jc(v,u)=c(v,u)+maxLY;%反向链路的流量裕度更新$ t- Z% n- c9 b# O: }: d* c
elseif flag(u,v)==1&&maxLY==c(u,v)%当这条边为前向边且是饱和边时 ' [, T+ j5 ?$ _ ?, aw(u,v)=BV;%更新权重值 l* C' a6 S% Q4 Vc(u,v)=c(u,v)-maxLY;%更新流量裕度值4 Z+ G1 r& j& G/ j+ u5 [
w(v,u)=-a(u,v);%反向链路权重更新' H# t3 T5 T5 j7 ]
elseif flag(u,v)==-1&&maxLY w(v,u)=a(v,u); * K; t b1 M/ Z" C6 rc(v,u)=c(v,u)+maxLY;3 x, m2 V: T/ L0 v# @7 C
w(u,v)=-a(v,u);, n% ~" M' t1 w& {. D/ G V
elseif flag(u,v)==-1&&maxLY==c(u,v)%当这条边为反向边且是饱和边时" D5 u% H/ ~$ F
w(v,u)=a(v,u);! o9 Y2 s, a/ f( C t
c(u,v)=c(u,v)-maxLY; 3 C1 f4 e9 ~; l6 Pw(u,v)=BV; ' K' @' s7 f2 B% U$ Jelse/ u5 n- }$ y0 [5 B
end+ f5 m3 u s$ w. c; Q
end 6 h. ] {# j, _$ E; o/ TMaxFlow2=sum(f(s,); ! y. D' P5 O/ r' vMinCost2=sum(sum(f.*a));" I2 N! n+ v% B2 l* g( d! z
if MaxFlow2<=V! ~3 e4 d" L t4 _; o" g- ?7 U
MaxFlow=MaxFlow2;1 j4 G( y+ ]' \; Z- ^. G: {' q# i0 c
MinCost=MinCost2;( d9 B4 Z( E; i. h6 h$ y
[L,R]=FLOYD(w,s,t); $ c9 |+ X _5 Z( D2 T' F) }/ [else 1 _/ k7 H3 q! R$ W! Lf=f1+prop*(f-f1);, Q) i' f" r3 o/ z- V$ h
MaxFlow=V;6 ?8 l4 n1 R& m/ z
MinCost=MinCost1+prop*(MinCost2-MinCost1);4 J$ W- H$ X" G3 y8 Z6 S4 T
return0 x. t8 M7 Z9 B5 T& Z# R. ?
end. U/ N6 ^# a+ h4 I; v4 m* B( u
if L(end) RE=1;%如果路径长度小于大数,说明路径存在8 z% b. P |5 ^& w& T
else ) l; a% i# C6 o) H& ERE=0;( P4 U. }- c$ y6 `( c; ~: Q5 i% @2 \
end 0 E- j! W( P8 ]4 P# eend . F7 D4 M$ u! D; {1 Mfunction [L,R]=FLOYD(w,s,t) / k1 n/ W5 c( V! o( cn=size(w,1);9 X( o& S) Z0 n. ^# M; Z
D=w;" s V+ c8 _' s" \, a
path=zeros(n,n);. y- B" ~2 r* e: D$ {% {* U8 X3 l
%以下是标准floyd算法4 K3 n. M4 L: I" l. i. E( Z. f
for i=1:n + f) m6 M% R, d1 F5 K& Vfor j=1:n ( T) I; a7 l5 c5 } g0 Q9 Y5 Vif D(i,j)~=inf " j" o- \( [+ `# d) f$ e/ opath(i,j)=j;+ J6 ]3 l: F9 I1 @' X$ I1 a
end ! \: v9 `6 Y" _6 x4 wend % X- H- z5 _: O2 J) i8 S4 @* u. J6 Vend & @' _: \. R9 e1 [% n& z, P* Cfor k=1:n; ?. M6 a3 d, _0 @# t
for i=1:n ) ]' F n6 @. h$ efor j=1:n 6 o6 J( J3 s0 [8 _ q6 S) q/ Qif D(i,k)+D(k,j) D(i,j)=D(i,k)+D(k,j);" F- J1 ^& N9 x$ I' N1 p$ R
path(i,j)=path(i,k); b6 n0 D0 H" ~' R
end9 F5 w% `0 c4 m6 p
end2 S% w3 _5 t$ V1 ~/ U6 A
end 2 C3 y/ a7 u, B+ Kend# s- D8 n7 s( B, }
L=zeros(0,0);% m7 c% K8 }2 B, P
R=s;$ A7 x9 m3 r; m4 B' e+ |
while 1 6 V; l# k( V0 f3 h' T6 ~, fif s==t7 p# p1 V+ C. \, d7 g# K3 I3 O t5 u
L=fliplr(L);9 L1 F+ {4 c& z0 l/ E' y5 `
L=[0,L];! p( {. m5 V+ U
return" @7 x1 k7 O# R' N( c- x1 Q
end : [. }) p8 e% e9 k* y5 }9 LL=[L,D(s,t)]; ! b: H. ]& s' F0 L* k1 m9 ?R=[R,path(s,t)];, h1 R0 Y( P) d# B
s=path(s,t);3 x$ A+ L% P- v, n0 }- L5 p- @
end