北京联合大学
实验报告
项目名称： 运筹学专题实验报告 学 院： 自动化 专 业： 物流工程 班 级： 1201B 学 号：[1**********]81 姓 名： 管水城 成 绩：
2015 年 5 月 6 日
实验三：使用matlab 求解最小费用最大流算问题
一、实验目的：
(1)使学生在程序设计方面得到进一步的训练; ，学习Matlab 语言进行程序设计求解最大流最小费用问题。
二、实验用仪器设备、器材或软件环境
计算机, Matlab R2006a
三、算法步骤、计算框图、计算程序等
1. 最小费用最大流问题的概念。
在网络D(V,A)中, 对应每条弧(vi,vj)IA,规定其容量限制为cij(cij\0),单位流量通过弧(vi,vj)的费用为dij(dij\0),求从发点到收点的最大流f, 使得流量的总费用d(f)为最小, 即mind(f)=E(vi,vj)IA
2. 求解原理。
若f 是流值为W 的所有可行流中费用最小者, 而P 是关于f 的所有可扩充链中费用最小的可扩充链, 沿P 以E 调整f 得到可行流fc, 则fc 是流值为(W+E)的可行流中的最小费用流。
根据这个结论, 如果已知f 是流值为W 的最小费用流, 则关键是要求出关于f 的最小费用的可扩充链. 为此, 需要在原网络D 的基础上构造一个新的赋权有向图E(f),使其顶点与D 的顶点相同, 且将D 中每条弧(vi,vj)均变成两个方向相反的弧(vi,vj)和(vj,vi)1新图E(f)中各弧的权值与f 中弧的权值有密切关系, 图E(f)中各弧的权值定义为：
新图E(f)中不考虑原网络D 中各个弧的容量cij. 为了使E(f)能比较清楚, 一般将长度为]的弧从图E(f)中略去. 由可扩充链费用的概念及图E(f)中权的定义可知, 在网络D 中寻求关于可行流f 的最小费用可扩充链, 等价于在图E(f)中寻求从发点到收点的最短路. 因图E(f)中有负权, 所以求E(f)中的最短路需用Floyd 算法。
1. 最小费用流算法的框图描述。
图一
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,k)+D(k,j);
R(i,j)=R(i,k);
end
end
end
k;
D;
R;
end
M=D(1,n);
3. 求解如下网络运输图中的最大流最小费用问题：
图2
打开matlab 软件，在COMND WINDOW窗口中输入矩阵程序如下：
n=5;
C=[0 10 8 0 0;0 0 0 2 7;0 5 0 10 0;0 0 0 0 4;0 0 0 0 0]
b=[0 4 1 0 0;0 0 0 6 1;0 2 0 3 0;0 0 0 0 2;0 0 0 0 0]
点击运行得到如下图：
图3
由上图实验结果可知，该问题的最大流为11，最小费用为55。
4. 求解如下最大流最小费用问题：
打开matlab 软件，在COMND WINDOW窗口中输入矩阵程序如下：
n=6;
C=[0 3 0 4 0 0;0 0 6 0 4 0;0 0 0 0 0 7;0 0 5 0 3 0;0 0 0 0 0 3;0 0 0 0 0 0]
b=[0 2 0 1 0 0;0 0 5 0 3 0;0 0 0 0 0 1;0 0 4 0 3 0;0 0 0 0 0 1;0 0 0 0 0 0]
点击运行得到如下图：
图4
由上图实验结果可知，该问题的最大流为7，最小费用为42。
四、实验总结
本实验在程序文件中所使用的计算最小费用最大流的算法并没有先用福德－富克逊法算出最大流，然后再用对偶法算出最小费用，而是将两种算法结合，最小费用和最大流一起算出。首先，福德－富克逊法要求对网络增加一个初始可行流，那么不妨设初始可行流为零流。然后再寻找增广链，可以采用对偶法以费用C 为权通过福德算法先找从起点至终点的最短路，再以该最短路为增广链调整流量，每一次调整都以矩阵a 记录调整的结果。为了能够满足增广链上正向弧非饱和、逆向弧非零流的条件，在每一次以C 为权寻找最短路之前，对费用C 矩阵进行调整。将正向饱和弧、逆向零流弧对应的C 值设为无穷大，非饱和弧的C 值设为初始值，这样一来，计算出的最短路径增广链就不会包括正向饱和弧与逆向零流弧了。每一次调整完网络流量之后，网络中的饱和弧、非饱和弧、零流弧会相互转化，因此要对网络中弧所对应的C 矩阵再次进行调整。调整的方法就是回到上边：将正向饱和弧、逆向零流弧对应的C 值设为无穷大、非饱和弧的C 值设为初始值，后再次以C 为权通过福德算法寻找最短路径，这样构成一个循环，直至以C 为权找不到一条从起点至终点的最短路径为止。找不到最短路径的标志就是福德算法返回从起点至终点的最短路径值为无穷大，此时网络已达最大流。