下面的最小费用最大流算法采用的是“基于Floyd最短路算法的Ford和Fulkerson迭加算法”,其基本思路为:把各条弧上单位流量的费用看成某种长度,用Floyd求最短路的方法确定一条自V1至Vn的最短路;再将这条最短路作为可扩充路,用求解最大流问题的方法将其上的流量增至最大可能值;而这条最短路上的流量增加后,其上各条弧的单位流量的费用要重新确定,如此多次迭代,最终得到最小费用最大流。) g4 ?' K# e! u2 {+ C F& X ]
O3 z! a( ]5 v* c' h- e+ H! L/ f& b; Kfunction [f,MinCost,MaxFlow]=MinimumCostFlow(a,c,V,s,t) 5 A3 O$ |/ F, z; Y% X6 C%% 基于Floyd最短路算法的Ford和Fulkerson迭加算法 1 T, c' x4 _# j7 [- ^%% 输入参数列表 4 u$ y# }; F, n: K* b1 T( A7 W& A R% a 单位流量的费用矩阵, ~, d# t8 C: J& N0 h
% c 链路容量矩阵2 I" I R4 ]6 ?/ m& s
% V 最大流的预设值,可为无穷大& i- ~0 B, @ v- Q( }$ L
% s 源节点 - w! F3 T7 O R) f( [, c% t 目的节点 3 ?, J9 p3 z$ Y6 [( c) k/ X' ^9 `$ w& m%% 输出参数列表+ U5 G1 p6 @! e9 b
% f 链路流量矩阵; p5 R3 B. u' Y
% MinCost 最小费用/ B& v0 z# |5 D( A' ]/ w
% MaxFlow 最大流量 ! Y9 b S2 f8 W%% 第一步:初始化2 V! ~/ E; j; U' b4 k9 s( p
N=size(a,1);%节点数目 + y3 s2 l& h7 s0 d ef=zeros(N,N);%流量矩阵,初始时为零流 ' I0 K: d: T* s: w2 @+ a$ q3 aMaxFlow=sum(f(s,);%最大流量,初始时也为零 % y) t2 x# m2 c6 W8 T' ]flag=zeros(N,N);%真实的前向边应该被记住- U1 w+ E5 T4 [
for i=1:N 3 V8 |# p' c% V: u! U( c/ Bfor j=1:N . q4 O6 ?1 n3 Y1 b. `/ |6 Vif i~=j&&c(i,j)~=04 d- j |+ X% Y
flag(i,j)=1;%前向边标记 " \% C' i# T( ]% Eflag(j,i)=-1;%反向边标记 . c8 b$ i6 d3 d! X! I' a4 e8 \end & d; T% O6 o, z; Z, p9 _) m- a9 ?if a(i,j)==inf$ V, K. h6 ?+ N7 c" \' }1 Z
a(i,j)=BV; % l! ]) Z! v8 v3 m+ w. B R: @w(i,j)=BV;%为提高程序的稳健性,以一个有限大数取代无穷大 7 h3 g# ]+ o* a7 z5 }) [end E) F s' E! L, l" j
end+ A# ~7 w4 s1 t B D0 R
end & [) V( |" d1 |5 w6 bif L(end) RE=1;%如果路径长度小于大数,说明路径存在0 L; E- ~8 A H. \
else ( u. u: R9 n) R/ D# VRE=0;$ S$ V9 V3 ^+ k; \5 j
end# D+ n8 U. c& a4 N1 b
%% 第二步:迭代过程6 U, {4 B+ p4 ~% u, z
while RE==1&&MaxFlow<=V%停止条件为达到最大流的预设值或者没有从s到t的最短路 , G3 [$ n9 C9 u7 b: {7 q%以下为更新网络结构# L( {( _$ x% Y$ k
MinCost1=sum(sum(f.*a)); * \9 f1 M0 L) h9 vMaxFlow1=sum(f(s,); 2 k" l4 a8 v' ]- Ef1=f;5 {0 y, v; n6 {5 K
TS=length(R)-1;%路径经过的跳数9 }" f+ Y0 ]# s+ h. G/ k
LY=zeros(1,TS);%流量裕度( ^; P- l, E7 P: d3 i; Y8 z
for i=1:TS$ u6 Q; U3 T0 p5 p
LY(i)=c(R(i),R(i+1)); ' b4 t5 F4 V. F) pend 0 K& c/ k" p o* pmaxLY=min(LY);%流量裕度的最小值,也即最大能够增加的流量 ) S' Z, Y1 ^& q) {) d5 sfor i=1:TS: f4 T' S* E; V$ Q6 n
u=R(i);9 o) T# o2 \& _7 s9 c
v=R(i+1); . T% a0 E) E$ w# ]2 C! ^if flag(u,v)==1&&maxLY f(u,v)=f(u,v)+maxLY;%记录流量值9 J/ S2 g/ ^" d
w(u,v)=a(u,v);%更新权重值( T5 K4 F: E' t+ K! d: I1 B
c(v,u)=c(v,u)+maxLY;%反向链路的流量裕度更新- t( w; c8 k( d- `
elseif flag(u,v)==1&&maxLY==c(u,v)%当这条边为前向边且是饱和边时 % R7 b6 m- Y7 ]; U4 y3 ww(u,v)=BV;%更新权重值9 A+ _" g3 q: D7 x9 k; [
c(u,v)=c(u,v)-maxLY;%更新流量裕度值, I0 F! U: b p7 i) M6 f
w(v,u)=-a(u,v);%反向链路权重更新5 n( r& e% H/ X; J% v
elseif flag(u,v)==-1&&maxLY w(v,u)=a(v,u);$ s2 H3 F# d0 t6 i# E/ Z
c(v,u)=c(v,u)+maxLY; $ m2 Y( |9 z6 \0 M8 }w(u,v)=-a(v,u);/ v, ^9 v) _6 H8 @' ]9 {2 L0 H
elseif flag(u,v)==-1&&maxLY==c(u,v)%当这条边为反向边且是饱和边时7 r$ X, C: c. V0 \, T
w(v,u)=a(v,u); 9 z' Z" l) A4 ]0 Y) W/ ic(u,v)=c(u,v)-maxLY;7 l4 z9 [+ F+ [! Z% L4 g7 C
w(u,v)=BV;, d6 b! W/ n. T$ V* M
else N3 I. ~7 Y- A5 `5 d& Zend 8 W+ z1 h8 v2 F7 I! G8 send 6 C. _$ L8 p; a9 yMaxFlow2=sum(f(s,);+ {/ ?% o/ G8 k c& Z" |0 R
MinCost2=sum(sum(f.*a));; I% j. e" Z: J5 ^
if MaxFlow2<=V . \2 V" v/ n8 H) Q6 N( CMaxFlow=MaxFlow2; ! {1 d4 m1 x* l: a; X5 dMinCost=MinCost2;1 u1 o3 I% N$ S% R, r: P
[L,R]=FLOYD(w,s,t); 1 `$ G( n0 s8 T- x0 Melse & Z" d! r. D( S: c. m6 Z/ p; mf=f1+prop*(f-f1); . A# _! @: [, ?* t7 {& lMaxFlow=V; # d( \7 |; |. x0 }! ]MinCost=MinCost1+prop*(MinCost2-MinCost1);" B/ a1 j6 j. p9 t# V% e
return 5 o% c) ^; ?( _# ^; c/ M. cend 1 {4 t: g; T/ b$ \+ m1 hif L(end) RE=1;%如果路径长度小于大数,说明路径存在( Q! J. I+ u, \
else 8 @6 ^7 Q% L. l- H+ Z/ g( ORE=0; 4 w8 e4 ?2 t5 I4 Lend + W1 t2 L& ]/ R9 `* g9 v6 \, uend % b) I, F5 Q1 w6 @. V. F0 W, Jfunction [L,R]=FLOYD(w,s,t) 4 K' S1 P9 f8 s) Q& p& J; B5 On=size(w,1); 5 h; H4 i3 g# N4 R4 GD=w;' Q# f4 M) Y% _$ V* p% F, Z; Y- F
path=zeros(n,n);8 B0 O2 G; A, B. V
%以下是标准floyd算法 2 w+ T1 b& W2 }7 efor i=1:n 3 w7 Y+ A- y2 C. x/ D8 x# M( gfor j=1:n( K: }+ P% @) T1 s& U8 I5 ?
if D(i,j)~=inf 3 O6 i( U# }5 T- k1 j, ~path(i,j)=j;0 q6 l0 U7 V3 t
end / G- T" |3 l2 W) send- S0 l# ^5 n2 c, W; y6 b; ~+ q
end 5 z8 p- c$ r3 [$ _' |' dfor k=1:n + @1 q! b( U+ J+ Y) u2 Q Dfor i=1:n8 |0 R! D# b, J( |: N0 H4 D
for j=1:n# A9 W5 S, E' @+ @
if D(i,k)+D(k,j) D(i,j)=D(i,k)+D(k,j); 7 M2 |* m# x( `6 ?3 d0 A) Kpath(i,j)=path(i,k);/ o3 z% B4 u/ W$ P+ _: K
end# C8 Y; b* l6 |9 B7 B
end ' S9 C% J) A% Kend( R- o E+ z6 b/ ^
end 6 g3 ~# P# h3 k) x/ m3 y4 |% D( ?L=zeros(0,0);% g# h/ u# n# y* U
R=s;1 h- w9 ^. `2 |3 ~
while 1 - L5 I( [% x8 z' V: Lif s==t. x, t- P/ c2 u2 s
L=fliplr(L); * r- _6 P: E8 x$ hL=[0,L]; 9 A Z1 d3 h1 ?$ J( w& p7 Z; }return c& N3 k* F/ z# \% jend5 z' ~, a6 I- m4 X$ \) p0 _. ]7 B
L=[L,D(s,t)]; , A. H* C, x" y0 z7 E( `5 C, qR=[R,path(s,t)]; 6 m2 b; }. v. J5 g2 { p( x, q+ Ls=path(s,t); # y% ^: T k! q& i, h& A, ], U+ xend