下面的最小费用最大流算法采用的是“基于Floyd最短路算法的Ford和Fulkerson迭加算法”,其基本思路为:把各条弧上单位流量的费用看成某种长度,用Floyd求最短路的方法确定一条自V1至Vn的最短路;再将这条最短路作为可扩充路,用求解最大流问题的方法将其上的流量增至最大可能值;而这条最短路上的流量增加后,其上各条弧的单位流量的费用要重新确定,如此多次迭代,最终得到最小费用最大流。6 o9 R4 Y! b' M) w4 H9 C7 a; r
" b" h7 l) r* h- Q* O- P) Ffunction [f,MinCost,MaxFlow]=MinimumCostFlow(a,c,V,s,t)& @0 x* B6 n. {; e$ B! D
%% 基于Floyd最短路算法的Ford和Fulkerson迭加算法9 w* V% t; F1 O% v! o6 p- S
%% 输入参数列表7 n# z" e y/ U# |% [
% a 单位流量的费用矩阵- D$ S) I" H8 a2 i( \* F
% c 链路容量矩阵. p: X( U; O1 o& z: m
% V 最大流的预设值,可为无穷大9 q! ?3 A; k4 n% P
% s 源节点$ ]. ^; p" s/ X6 E
% t 目的节点 / U% g( K8 i3 P4 S2 `- T: E3 u%% 输出参数列表 ) q4 r: O' k) E% U4 m% f 链路流量矩阵 8 A" y1 u/ C& d1 ^) C% MinCost 最小费用2 r+ b G; \' m( ^ y, G$ a
% MaxFlow 最大流量' Q4 U7 T- k" n
%% 第一步:初始化 + M z- \0 K; {9 tN=size(a,1);%节点数目 4 X' E& t) x) y9 b5 of=zeros(N,N);%流量矩阵,初始时为零流. P# T K# Z0 J4 d
MaxFlow=sum(f(s,);%最大流量,初始时也为零0 [% `; [2 ]' G" y$ @7 f
flag=zeros(N,N);%真实的前向边应该被记住2 t! m* @2 b( c# S. z: R
for i=1:N 6 h# m! o3 t, [for j=1:N 3 ^% K. Z9 x; m% W9 Aif i~=j&&c(i,j)~=0 3 K: H& S+ B+ |/ q: Rflag(i,j)=1;%前向边标记 ( v# M& g4 h- N9 _- _7 K; P% S: a0 rflag(j,i)=-1;%反向边标记' a; S0 p# F9 x
end . l E! g0 g/ n8 L& K! zif a(i,j)==inf 9 M0 X# v# G& A" p. W" ka(i,j)=BV; 1 _! Q# g7 `5 y( v* y+ Tw(i,j)=BV;%为提高程序的稳健性,以一个有限大数取代无穷大 " d5 i: T4 \- [; \end " @ ^" M6 Z( [; _4 L0 hend7 _1 D7 f! U- J8 i# O2 d
end$ I" _7 @ g: m- `0 r! a
if L(end) RE=1;%如果路径长度小于大数,说明路径存在 % o3 t9 z1 H. g- Kelse % U: [' A6 m, I# xRE=0; 9 t+ u+ L7 _: Send ) Z- I9 Y" m6 ^%% 第二步:迭代过程 # u" \/ I) O6 N# l. X. f" Iwhile RE==1&&MaxFlow<=V%停止条件为达到最大流的预设值或者没有从s到t的最短路) C2 M& @! z' l" W/ Q6 g7 ?- _
%以下为更新网络结构4 C- }. M7 n# W$ F
MinCost1=sum(sum(f.*a)); : x8 X: u5 D0 X* E7 [MaxFlow1=sum(f(s,);* ^3 U. Q8 k7 R& x4 W8 u! S! }
f1=f; & N9 p2 M! Q5 y' hTS=length(R)-1;%路径经过的跳数 ( U P- ^' E# m9 b% H4 v1 NLY=zeros(1,TS);%流量裕度 9 y. A% l0 n: W" z, H: \. ufor i=1:TS . P% T+ ?7 `% i ?! h! nLY(i)=c(R(i),R(i+1));% W; ?. m6 [" ^) T
end! I+ a( h4 J% }- O
maxLY=min(LY);%流量裕度的最小值,也即最大能够增加的流量& o1 K. \. M8 h1 q; [) ?
for i=1:TS % C# ~3 {/ j" s! `% T. ou=R(i);3 F( j! h$ x( T+ u9 i
v=R(i+1);* V, [2 ?: m K/ R6 U+ D
if flag(u,v)==1&&maxLY f(u,v)=f(u,v)+maxLY;%记录流量值 ( v, `6 {; @0 Tw(u,v)=a(u,v);%更新权重值 * I1 z$ Y7 i3 S3 dc(v,u)=c(v,u)+maxLY;%反向链路的流量裕度更新, p% G+ h& [- G
elseif flag(u,v)==1&&maxLY==c(u,v)%当这条边为前向边且是饱和边时) \' f; ` [' p0 I5 K
w(u,v)=BV;%更新权重值 M( O- S# t# z1 k4 Pc(u,v)=c(u,v)-maxLY;%更新流量裕度值% A( o6 `' P( y0 J! E P
w(v,u)=-a(u,v);%反向链路权重更新 % l$ L! B+ h7 J( {4 uelseif flag(u,v)==-1&&maxLY w(v,u)=a(v,u); , j1 v3 l" ?$ ic(v,u)=c(v,u)+maxLY;7 M2 T# V G7 v/ R
w(u,v)=-a(v,u); 0 M0 s a: j! {elseif flag(u,v)==-1&&maxLY==c(u,v)%当这条边为反向边且是饱和边时5 F# U* O) |6 x8 S3 S/ _% W
w(v,u)=a(v,u);8 E' ?$ d4 p# d8 } c$ Y+ N# {7 T I
c(u,v)=c(u,v)-maxLY; 8 K& b) B4 }2 K! y- r" Yw(u,v)=BV;) e. O7 t9 n5 L8 }
else # J7 G2 p& T7 E: K) vend - z$ @ R9 I; A4 Iend 2 M0 U; m' ?8 D2 WMaxFlow2=sum(f(s,); 1 b1 b" S0 s% u: ~! W" k+ d# uMinCost2=sum(sum(f.*a)); 0 f3 j+ i$ g7 \6 ~if MaxFlow2<=V * ]4 w9 s- ?9 Y1 b8 O" ]MaxFlow=MaxFlow2; # v. N7 V4 o; j& ~9 J. OMinCost=MinCost2; & f5 |( W; {, R8 U[L,R]=FLOYD(w,s,t); ! p8 M7 X) G% P2 r6 ielse p/ ` W/ i5 R: V- [! k( O1 z, C
f=f1+prop*(f-f1);: P. ?! g! S2 P% G w3 x/ k, t- Y
MaxFlow=V; " d# k* F5 N( l5 ]* h) pMinCost=MinCost1+prop*(MinCost2-MinCost1); 1 I! F8 Z5 p( mreturn$ C3 K6 S3 v# d, C" p
end" h3 H: O V3 B( [, D( F' y' a) b
if L(end) RE=1;%如果路径长度小于大数,说明路径存在9 p2 j, Q; i% K$ i3 \
else 3 t( R% d2 K* S2 x$ qRE=0;2 H; u) B! C9 U* v6 d* ~
end8 L) ^( J, l( ?1 m8 a9 P+ _) [
end2 I) N) g9 k. I9 o: I( f$ y0 `
function [L,R]=FLOYD(w,s,t) ' t- r* X/ B1 W! fn=size(w,1); - Z4 [# \! o# J( Y6 lD=w;# z* r1 [& @2 E, D# |0 y% Y3 q
path=zeros(n,n); 3 r3 R% S0 g* `%以下是标准floyd算法 {6 X8 x1 s" `6 B0 M( |) Z- \4 g
for i=1:n ; t1 y9 ^+ N1 c! O) B. R( efor j=1:n7 ^ v, h2 D. {
if D(i,j)~=inf 0 N* `' e# i4 a3 k# r$ kpath(i,j)=j; 7 I* T/ h( B$ V1 D; o) Aend . D0 G. l+ M+ N, Nend 8 O# c9 [, o3 jend " [% a$ O, t9 h* X3 r5 F6 n: r8 h2 Efor k=1:n , N8 ~ }) a" h8 K6 _for i=1:n% N$ ]/ k7 X: M; R( X
for j=1:n ) n% b: w/ w/ yif D(i,k)+D(k,j) D(i,j)=D(i,k)+D(k,j);! D6 b- P6 k3 h7 f, W6 ?; l+ e
path(i,j)=path(i,k); : f- z: K3 H& v* c( J- xend 4 u- P8 B/ o' d& ^2 c' C! [0 V+ b# xend/ t4 R$ G/ H/ l( V+ T, \# B: Z+ F
end ; S3 @) w' ?( O3 H. Fend1 o' Z- I8 z* z2 H1 m
L=zeros(0,0);8 r" d5 k' X+ u$ N* Y3 _
R=s;# L! w& ~* M0 {3 L
while 1 / g- r: E% h d9 O0 \; c2 Lif s==t * d/ r0 P5 A+ F/ o f/ KL=fliplr(L);5 O4 Q0 X( y% K6 d8 l( `3 @
L=[0,L];0 h3 b2 V) E0 b3 Q/ x
return+ v( Z& b. c1 U4 b. A3 Y, L
end# Y5 a- ]' k5 g( |* b! ]. v; c* K$ X! H
L=[L,D(s,t)];/ X" Y* g: Y) u# w$ X/ A" H
R=[R,path(s,t)]; ; c* Q0 k: g# e2 Z7 Fs=path(s,t);9 Q9 m# p, x7 L S
end