它山之拾

我在程序文件中所使用的计算最小费用最大流的算法并没有先用福德－富克逊法算出最大流，然后再用对偶法算出最小费用，而是将两种算法结合，最小费用和最大流一起算出。
首先，福德－富克逊法要求对网络增加一个初始可行流，那么不妨设初始可行流为零流。然后再寻找增广链，可以采用对偶法以费用C为权通过福德算法先找从起点至终点的最短路，再以该最短路为增广链调整流量，每一次调整都以矩阵a记录调整的结果。为了能够满足增广链上正向弧非饱合、逆向弧非零流的条件，在每一次以C为权寻找最短路之前，对费用C矩阵进行调整。将正向饱合弧、逆向零流弧对应的C值设为无穷大，非饱合弧的C值设为初始值，这样一来，计算出的最短路径增广链就不会包括正向饱合弧与逆向零流弧了。每一次调整完网络流量之后，网络中的饱合弧、非饱合弧、零流弧会相互转化，因此要对网络中弧所对应的C矩阵再次进行调整。调整的方法就是回到上边：将正向饱合弧、逆向零流弧对应的C值设为无穷大、非饱合弧的C值设为初始值，后再次以C为权通过福德算法寻找最短路径，这样构成一个循环，直至以C为权找不到一条从起点至终点的最短路径为止。找不到最短路径的标志就是福德算法返回从起点至终点的最短路径值为无穷大，此时网络已达最大流。可以见下面：





使用上面的算法，有关最大流可能会提出这样两个问题：
1.     由于是以C为权值寻找的最短路径，那么是不是不能穷尽从起点至终点的路径，如果不能穷尽那就有可能存在没有找到的增广链，还能找到增广链那就不是最大流？
由于这个算法在寻找最短路径前对C矩阵进行了调整，将非增广链的路径排除在寻找范围之外，所以它保证找到的每一条路径都是增广链。而且它每次得到的都是剩余增广链中的最短增广链，而且找到一条就令其花费为无穷大，以至下次寻找到的是次短增广链。而增广链的数量又是有限的，这样一直寻找下去就会将所有增广链穷尽。就好比有1~10十个数字，每次取剩余中的最小数字，取后不放回，最终会将10个数字全部取完。既然增广链可以穷尽，当找不到增广链时，网络流量即达最大流量，这是福德－富克逊法的思想。穷尽之后C矩阵中起点至终点的各条路径值皆为无穷大，因此福德算法会返回从起点到终点的最短路径值为无穷大，可以将此做为退出循环的标志。
2.     会不会已经达到了网络最大流，但该程序还继续循环，导致错误？
程序已经将不符合增广链要求的路径排除在寻找范围之外，于是程序继续循环证明还能找到增加广链，能找到增广链就没有达到最大流。
对于最小费用问题，由于该程序在计算最小费用方面的算法是优先选择费用低的链路，当费用低的链路在保证网络最大流的前提下不能再提供流量时，再转到次低费用的链路上，与“对偶法”的思想一致。因此，最小费用的计算也是正确的。

计算最小费用最大流MATLAB源代码,文件名为mp_mc.m
function[Mm,mc,Mmr]=mp_mc(a,c)
A=a; %各路径最大承载流量矩阵
C=c; %各路径花费矩阵
Mm=0; %初始可行流设为零
mc=0; %最小花费变量
mcr=0;
mrd=0;
n=0;
while mrd~=inf %一直叠代到以花费为权值找不到最短路径
     for i=1:(size(mcr',1)-1)
           if a(mcr(i),mcr(i+1))==inf
                 ta=A(mcr(i+1),mcr(i))-a(mcr(i+1),mcr(i));
           else
                 ta=a(mcr(i),mcr(i+1));
           end
           n=min(ta,n); %将最短路径上的最小允许流量提取出来
     end
     for i=1:(size(mcr',1)-1)
           if a(mcr(i),mcr(i+1))==inf
                 a(mcr(i+1),mcr(i))=a(mcr(i+1),mcr(i))+n;
           else
                 a(mcr(i),mcr(i+1))=a(mcr(i),mcr(i+1))-n;
           end
     end
     Mm=Mm+n; %将每次叠代后增加的流量累加，叠代完成时就得到最大流量
     for i=1:size(a,1)
           for j=1:size(a',1)
                 if i~=j&a(i,j)~=inf
                       if a(i,j)==A(i,j) %零流弧
                             c(j,i)=inf;
                             c(i,j)=C(i,j);
                       elseif a(i,j)==0 %饱合弧
                             c(i,j)=inf;
                             c(j,i)=C(j,i);
                       elseif a(i,j)~=0 %非饱合弧
                             c(j,i)=C(j,i);
                             c(i,j)=C(i,j);
                       end
                 end
           end
     end
     [mcr,mrd]=floyd_mr(c) %进行叠代，得到以花费为权值的最短路径矩阵(mcr)和数值(mrd)
     n=inf;
end
%下面是计算最小花费的数值
for i=1:size(A,1)
     for j=1:size(A',1)
           if A(i,j)==inf
                 A(i,j)=0;
           end
           if a(i,j)==inf
                 a(i,j)=0;
           end
     end
end
Mmr=A-a; %将剩余空闲的流量减掉就得到了路径上的实际流量，行列交点处的非零数值就是两点间路径的实际流量
for i=1:size(Mmr,1)
     for j=1:size(Mmr',1)
           if Mmr(i,j)~=0
                 mc=mc+Mmr(i,j)*C(i,j); %最小花费为累加各条路径实际流量与其单位流量花费的乘积
           end
     end
end


利用福得算法计算最短路径MATLAB源代码,文件名为floyd_mr.m
function[mr,mrd]=floyd_mr(a)
n=size(a,1);
[D,R]=floyd(a); %通过福德算法得到距离矩阵(D)和路径矩阵(R)
mrd=D(1,n); %提取从起点1到终点n的最短距离
rd=R(1,n); %提取从起点1开始沿最短路径上下一个点的编号(rd)
mr=[1,rd]; %从起点1开始沿最短路径到rd点的最短路径

while rd~=n %通过循环将最短路径依次提取出来，直到rd点就是最后一个点
     mr=[mr,R(rd,n)];
     rd=R(rd,n);
end

福得算法MATLAB源代码,文件名为floyd.m
function[D,R]=floyd(a)

n=size(a,1);
D=a;
for i=1:n
  for j=1:n
  R(i,j)=j;
  end
end
R;

for k=1:n
  for i=1:n
    for j=1:n
      if D(i,k)+D(k,j)<D(i,j)
        D(i,j)=D(i,k)+D(k,j);
        R(i,j)=R(i,k);
      end
    end
  end
  k;
  D;
  R;
end
M=D(1,n);