图论算法及其MATLAB程序代码.pdf

图论算法及其 MATLAB程序代码 求赋权图 G=(V,E,F)中任意两点间的最短路的 Warshall-Floyd 算法：设 A=(aij)nn为赋权图 G=(V,E,F)的矩阵,当 vivjE 时 aij=F(vivj),否则取 aii=0,aij=+(ij),dij表示从 vi到 vj点的距离,rij表示从 vi到 vj点的最短路中一个点的编号.赋初值.对所有 i,j,dij=aij,rij=j.k=1.转向 更新 dij,rij.对所有 i,j,若 dik+dk jdij,则令 dij=dik+dk j,rij=k,转向.终止判断.若 dii0,则存在一条含有顶点 vi的负回路,终止;或者 k=n 终止;否则令 k=k+1,转向.最短路线可由 rij得到.例 1 求图 6-4 中任意两点间的最短路.解：用 Warshall-Floyd 算法,MATLAB程序代码如下：n=8;A=0 2 8 1 Inf Inf Inf Inf 2 0 6 Inf 1 Inf Inf Inf 8 6 0 7 5 1 2 Inf 1 Inf 7 0 Inf Inf 9 Inf Inf 1 5 Inf 0 3 Inf 8 Inf Inf 1 Inf 3 0 4 6 Inf Inf 2 9 Inf 4 0 3 Inf Inf Inf Inf 8 6 3 0;%MATLAB 中,Inf表示 D=A;%赋初值 for(i=1:n)for(j=1:n)R(i,j)=j;end;end%赋路径初值 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);%更新 dij R(i,j)=k;end;end;end%更新 rij k%显示迭代步数 D%显示每步迭代后的路长 R%显示每步迭代后的路径 pd=0;for i=1:n%含有负权时 if(D(i,i)0)x(k)=A(i,j);%数组 x记录 A中不同的正数 kk=1;%临时变量 for(s=1:k-1)if(x(k)=x(s)kk=0;break;end;end%排除相同的正数 k=k+kk;end;end;end k=k-1%显示 A中所有不同正数的个数 for(i=1:k-1)for(j=i+1:k)%将 x中不同的正数从小到大排序 if(x(j)0)kk=kk+1;zz=z;end;end%寻找 TT 中的树枝 if(kk=1)TT(y,zz)=0;TT(zz,y)=0;pd=0;end;end%砍掉 TT 中的树枝 if(pd)break;end;end%已砍掉了 TT 中所有的树枝 pd=0;%判断 TT 中是否有圈 for(y=1:n-1)for(z=y+1:n)if(TT(y,z)0)pd=1;break;end;end;end if(pd)T(i,j)=0;T(j,i)=0;%假如 TT 中有圈 else q=q+1;end;end;end;end;end T%显示近似最小生成树 T,程序结束 求二部图 G 的最大匹配的算法(匈牙利算法),其基本思想是：从 G 的任意匹配 M 开始,对 X 中所有 M 的非饱和点,寻找 M 增广路.若不存在 M 增广路,则 M 为最大匹配;若存在 M 增广路 P,则将 P 中 M 与非 M 的边互换得到比 M 多一边的匹配 M1,再对 M1重复上述过程.设 G=(X,Y,E)为二部图,其中 X=x1,x2,xn,Y=y1,y2,yn.任取 G 的一初始匹配 M(如任取 eE,则 M=e是一个匹配).令 S=,T=,转向.若 M 饱和 X S 的所有点,则 M 是二部图 G 的最大匹配.否则,任取 M 的非饱和点uX S,令 S=S u,转向.记 N(S)=v|uS,uvE.若 N(S)=T,转向.否则取 yN(S)T.若 y是 M的饱和点,转向,否则转向.设 x yM,则令 S=S x,T=T y,转向.u y 路是 M增广路,设为 P,并令 M=MP,转向.这里 MP=MP MP,是对称差.由于计算 M增广路 P 比较麻烦,因此将迭代步骤改为：将 X 中 M 的所有非饱和点(不是 M 中某条边的端点)都给以标号 0 和标记*,转向.若 X 中所有有标号的点都已去掉了标记*,则 M 是 G 的最大匹配.否则任取 X 中一个既有标号又有标记*的点 xi,去掉 xi的标记*,转向.找出在 G 中所有与 xi邻接的点 yj(即 xi yjE),若所有这样的 yj都已有标号,则转向,否则转向.对与xi邻接且尚未给标号的yj都给定标号i.若所有的yj都是M的饱和点,则转向,否则逆向返回.即由其中M的任一个非饱和点yj的标号i找到xi,再由xi的标号k找到yk,最后由 yt的标号 s找到标号为 0 的 xs时结束,获得M 增广路 xs yt xi yj,记P=xs yt,xi yj,重新记 M 为 MP,转向.将 yj在 M 中与之邻接的点 xk(即 xk yjM),给以标号 j 和标记*,转向.例 1 求图 6-9 中所示的二部图 G 的最大匹配.匈牙利算法的 MATLAB程序代码如下：m=5;n=5;A=0 1 1 0 0 1 1 0 1 1 0 1 1 0 0 0 1 1 0 0 0 0 0 1 1;M(m,n)=0;for(i=1:m)for(j=1:n)if(A(i,j)M(i,j)=1;break;end;end%求初始匹配 M if(M(i,j)break;end;end%获得仅含一条边的初始匹配 M while(1)for(i=1:m)x(i)=0;end%将记录 X中点的标号和标记*for(i=1:n)y(i)=0;end%将记录 Y中点的标号和标记*for(i=1:m)pd=1;%寻找 X中 M的所有非饱和点 for(j=1:n)if(M(i,j)pd=0;end;end if(pd)x(i)=-n-1;end;end%将 X中 M的所有非饱和点都给以标号 0 和标记*,程序中用 n+1 表示 0 标号,标号为负数时表示标记*pd=0;while(1)xi=0;for(i=1:m)if(x(i)1)k=k-1;for(j=1:k)pdd=1;for(i=1:m)if(M(i,yy(j)x(i)=-yy(j);pdd=0;break;end;end%将yj在 M中与之邻接的点 xk(即 xkyjM),给以标号 j 和标记*if(pdd)break;end;end if(pdd)k=1;j=yy(j);%yj 不是 M的饱和点 while(1)P(k,2)=j;P(k,1)=y(j);j=abs(x(y(j);%任取 M的一个非饱和点 yj,逆向返回 if(j=n+1)break;end%找到 X中标号为 0 的点时结束,获得 M-增广路 P k=k+1;end for(i=1:k)if(M(P(i,1),P(i,2)M(P(i,1),P(i,2)=0;%将匹配 M 在增广路 P 中出现的边去掉 else M(P(i,1),P(i,2)=1;end;end%将增广路 P 中没有在匹配 M中出现的边加入到匹配 M中 break;end;end;end if(pd)break;end;end%假如 X中所有有标号的点都已去掉了标记*,算法终止 M%显示最大匹配 M,程序结束 图 6-9 利用可行点标记求最佳匹配的算法步骤如下：设 G=(X,Y,E,F)为完备的二部赋权图,L 是其一个初始可行点标记,通常取.,0)(,|)(max)(YyXxyLYyxyFxL=M 是 GL的一个匹配.若X 的每个点都是 M 的饱和点,则 M 是最佳匹配.否则取 M 的非饱和点 uX,令 S=u,T=,转向.记NL(S)=v|uS,uvEL.若NL(S)=T,则GL没有完美匹配,转向.否则转向.调整可行点标记,计算 aL=min L(x)+L(y)F(x y)|xS,yY T.由此得新的可行顶点标记 H(v)=,),(,)(,)(TvSvvLavLavLLL+令 L=H,GL=GH,重新给出 GL的一个匹配 M,转向.取 yNL(S)T,若 y 是 M 的饱和点,转向.否则,转向.设 x yM,则令 S=S x,T=T y,转向.在 GL中的 u y 路是 M 增广路,记为 P,并令 M=MP,转向.利用可行点标记求最佳匹配算法的 MATLAB程序代码如下：n=4;A=4 5 5 1 2 2 4 6 4 2 3 3 5 0 2 1;for(i=1:n)L(i,1)=0;L(i,2)=0;end for(i=1:n)for(j=1:n)if(L(i,1)L(S(i),1)+L(j,2)-A(S(i),j)al=L(S(i),1)+L(j,2)-A(S(i),j);end;end;end for(i=1:jss)L(S(i),1)=L(S(i),1)-al;end%调整可行点标记 for(j=1:jst)L(T(j),2)=L(T(j),2)+al;end%调整可行点标记 for(i=1:n)for(j=1:n)%生成子图 GL if(L(i,1)+L(j,2)=A(i,j)Gl(i,j)=1;else Gl(i,j)=0;end M(i,j)=0;k=0;end;end ii=0;jj=0;for(i=1:n)for(j=1:n)if(Gl(i,j)ii=i;jj=j;break;end;end if(ii)break;end;end%获得仅含 Gl 的一条边的初始匹配 M M(ii,jj)=1;break else%NL(S)T for(j=1:jsn)pd=1;%取 yNL(S)T for(k=1:jst)if(T(k)=NlS(j)pd=0;break;end;end if(pd)jj=j;break;end;end pd=0;%判断 y是否为 M的饱和点 for(i=1:n)if(M(i,NlS(jj)pd=1;ii=i;break;end;end if(pd)jss=jss+1;S(jss)=ii;jst=jst+1;T(jst)=NlS(jj);%S=Sx,T=Ty else%获得 Gl 的一条 M-增广路,调整匹配 M for(k=1:jst)M(S(k),T(k)=1;M(S(k+1),T(k)=0;end if(jst=0)k=0;end M(S(k+1),NlS(jj)=1;break;end;end;end;end MaxZjpp=0;for(i=1:n)for(j=1:n)if(M(i,j)MaxZjpp=MaxZjpp+A(i,j);end;end;end M%显示最佳匹配 M MaxZjpp%显示最佳匹配 M的权,程序结束 从一个可行流 f 开始,求最大流的 Ford-Fulkerson 标号算法的基本步骤：标号过程 给发点 vs以标号(+,+),s=+.选择一个已标号的点 x,对于 x 的所有未给标号的邻接点 y,按下列规则处理：当 yxE,且 f yx 0 时,令 y=min f yx,x,并给 y 以标号(x ,y).当 xyE,且 f xyC xy时,令 y=min C xy f xy,x,并给 y 以标号(x+,y).重复直到收点 vt被标号或不再有点可标号时为止.若 vt得到标号,说明存在一条可增广链,转调整过程;若 vt未得到标号,标号过程已无法进行时,说明 f 已经是最大流.调整过程 决定调整量=vt,令 u=vt.若 u 点标号为(v+,u),则以f vu+代替f vu;若 u 点标号为(v,u),则以 f vu 代替 f vu.若 v=vs,则去掉所有标号转重新标号;否则令 u=v,转.算法终止后,令已有标号的点集为 S,则割集(S,S c)为最小割,从而 Wf=C(S,S c).例 1 求图 6-19 所示网络的最大流.利用 Ford-Fulkerson 标号法求最大流算法的 MATLAB 程序代码如下：n=8;C=0 5 4 3 0 0 0 0 0 0 0 0 5 3 0 0 0 0 0 0 0 3 2 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0;%弧容量 for(i=1:n)for(j=1:n)f(i,j)=0;end;end%取初始可行流 f 为零流 for(i=1:n)No(i)=0;d(i)=0;end%No,d 记录标号 图 6-19 while(1)No(1)=n+1;d(1)=Inf;%给发点 vs 标号 while(1)pd=1;%标号过程 for(i=1:n)if(No(i)%选择一个已标号的点 vi for(j=1:n)if(No(j)=0&f(i,j)d(i)d(j)=d(i);end elseif(No(j)=0&f(j,i)0)%对于未给标号的点 vj,当 vjvi 为非零流弧时 No(j)=-i;d(j)=f(j,i);pd=0;if(d(j)d(i)d(j)=d(i);end;end;end;end;end if(No(n)|pd)break;end;end%若收点 vt 得到标号或者无法标号,终止标号过程 if(pd)break;end%vt 未得到标号,f 已是最大流,算法终止 dvt=d(n);t=n;%进入调整过程,dvt 表示调整量 while(1)if(No(t)0)f(No(t),t)=f(No(t),t)+dvt;%前向弧调整 elseif(No(t)0&f(i,j)=0)a(i,j)=b(i,j);elseif(C(i,j)0&f(i,j)=C(i,j)a(j,i)=-b(i,j);elseif(C(i,j)0)a(i,j)=b(i,j);a(j,i)=-b(i,j);end;end;end for(i=2:n)p(i)=Inf;s(i)=i;end%用 Ford 算法求最短路,赋初值 for(k=1:n)pd=1;%求有向赋权图中 vs 到 vt 的最短路 for(i=2:n)for(j=1:n)if(p(i)p(j)+a(j,i)p(i)=p(j)+a(j,i);s(i)=j;pd=0;end;end;end if(pd)break;end;end%求最短路的 Ford 算法结束 if(p(n)=Inf)break;end%不存在 vs 到 vt 的最短路,算法终止.注意在求最小费用最大流时构造有向赋权图中不会含负权回路,所以不会出现 k=n dvt=Inf;t=n;%进入调整过程,dvt 表示调整量 while(1)%计算调整量 if(a(s(t),t)0)dvtt=C(s(t),t)-f(s(t),t);%前向弧调整量 elseif(a(s(t),t)dvtt)dvt=dvtt;end if(s(t)=1)break;end%当 t 的标号为 vs 时,终止计算调整量 t=s(t);end%继续调整前一段弧上的流 f pd=0;if(wf+dvt=wf0)dvt=wf0-wf;pd=1;end%如果最大流量大于或等于预定的流量值 t=n;while(1)%调整过程 if(a(s(t),t)0)f(s(t),t)=f(s(t),t)+dvt;%前向弧调整 elseif(a(s(t),t)0)f(t,s(t)=f(t,s(t)-dvt;end%后向弧调整 if(s(t)=1)break;end%当 t 的标号为 vs 时,终止调整过程 t=s(t);end if(pd)break;end%如果最大流量达到预定的流量值 wf=0;for(j=1:n)wf=wf+f(1,j);end;end%计算最大流量 zwf=0;for(i=1:n)for(j=1:n)zwf=zwf+b(i,j)*f(i,j);end;end%计算最小费用 f%显示最小费用最大流 图 6-22 wf%显示最小费用最大流量 zwf%显示最小费用,程序结束