G = graph(s,t,weights) ——使用节点对组s,t和权重向量weights
创建赋权无向图
G = graph(s,t,weights,nodes)——使用字符向量元胞数组nodes
指定节点名称
G = graph(s,t,weights,num) ——使用数值标量num指定图中的
节点数
G = graph(A[,nodes],type)——仅使用A的上或下三角形阵构造
赋权图,type可以是’upper’ 或
‘lower’
sparse邻接矩阵转系数矩阵
例1无权无向图
clc, clear, close all
E=[1,2;1,3;2,3;3,2;3,5;4,2;4,6;5,2;5,4;6,5]
s = E(:,1); t = E(:,2);
nodes =cellstr(strcat('v',int2str([1:6]')))
G = digraph(s, t, [], nodes);
plot(G)
W1 = adjacency(G) %邻接矩阵的系数矩阵
W2 = incidence(G)%关联矩阵的系数矩阵
其中clear all:清除工作空间的所有变量,函数,和MEX文件,strcat
cellstr
邻接矩阵表示结点与结点之间的关系,关联矩阵表示结点与边的关系
例2有权无向图最短路
clc, clear, close all
E = [1,2,50; 1,4,40; 1,5,25; 1,6,10; 2,3,15; 2,4,20; 2,6,25;
3,4,10;3,5,20; 4,5,10; 4,6,25; 5,6,55];
G = graph(E(:,1), E(:,2), E(:,3));
p = plot(G,'EdgeLabel',G.Edges.Weight,'Layout', 'circle');
[path2, d2] = shortestpath(G, 1, 2) %求1到2的最短路
highlight(p,path2,'EdgeColor','r','LineWidth',1.5)%在图中画出path2的路线
[path3, d3] = shortestpath(G, 1, 3, 'method','positive') %同上,'method','positive'可省略
highlight(p,path3,'EdgeColor', 'm','LineWidth',1.5)
另一个例子
clc, clear, close all
a(1,2)=2;a(1,3)=8;a(1,4)=1;a(2,3)=6;a(2,5)=1;
a(3,4)=7;a(3,5)=5;a(3,6)=1;a(3,7)=2;a(4,7)=9;
a(5,6)=3;a(5,8)=2;a(5,9)=9;a(6,7)=4;a(6,9)=6;
a(7,9)=3;a(7,10)=1;a(8,9)=7;a(8,11)=9;
a(9,10)=1;a(9,11)=2;a(10,11)=4;
a=a';
[i,j,v]=find(a);
b=sparse(i,j,v,11,11);
[x,y,z]=graphshortestpath(b,1,11,'directed',false)
其中z表示1(graphshorestpath里设定的,不是默认)到每个点最短路的倒数第二个节点,比如z3=6,1到3的最短路从1到2到5到6再到3
从两个例子可以看出shortestpath对已经生成的图使用,graphshortestpath对稀疏矩阵使用,下面例7的graphmaxflow和maxflow同理
例3 Dijkstra
clc, clear, a=zeros(6);
a(1,[2,4,5,6])=[50,40,25,10];
a(2,[3,4,6])=[15,20,25];
a(3,[4,5])=[10,20];
a(4,[5,6])=[10,25];
a(5,6)=55;
G = graph(a,'upper');
d = distances(G,'Method','positive')
例4 Prime求最小生成树(好像用不到)
a=zeros(7);
a(1,2)=50; a(1,3)=60; a(2,4)=65; a(2,5)=40; a(3,4)=52;
a(3,7)=45; a(4,5)=50; a(4,6)=30;a(4,7)=42; a(5,6)=70;
a=a+a';
a(find(a==0))=inf;
result=[];p=1;tb=2:length(a);
while size(result,2)~=length(a)-1
temp=a(p,tb);
temp=temp(:);
d=min(temp);
[jb,kb]=find(a(p,tb)==d,1);
j=p(jb);k=tb(kb);
result=[result,[j;k;d]];p=[p,k];tb(find(tb==k))=[];
end
result
例5 Kruskal求最小生成树(好像用不到)
a(1,2)=50; a(1,3)=60; a(2,4)=65; a(2,5)=40;a(3,4)=52;a(3,7)=45;
a(4,5)=50; a(4,6)=30;a(4,7)=42; a(5,6)=70;
[i,j,b]=find(a);data=[i';j';b'];index=data(1:2,:);
loop=length(a)-1;result=[];
while length(result)<loop
temp=min(data(3,:));
flag=find(data(3,:)==temp);
flag=flag(1);
v1=index(1,flag);
v2=index(2,flag);
if v1~=v2
result=[result,data(:,flag)];
end
index(find(index==v2))=v1;
data(:,flag)=[];index(:,flag)=[];
end
result
例6最小生成树函数
clc, clear, close all, a=zeros(7);
a(1,[2 3])=[50,60];a(3,[4,7])=[52,45];
a(2,[4 5])=[65 40]; a(4,[5:7])=[50,30,42];
a(5,6)=70;
s=cellstr(strcat('v',int2str([1:7]')));
G=graph(a,s,'upper');
p=plot(G,'EdgeLabel',G.Edges.Weight)
T=minspantree(G,'Method','sparse')
L=sum(T.Edges.Weight), highlight(p,T)
其中sum为最小生成树的权重,常用于连通每个点,求最小权值,upper表示输入的矩阵为上三角,反之lower
另一个例子
clc, clear, xy=load('data4_12.txt'); d=mandist(xy); %求xy的两两列向量间的绝对值距离
s=cellstr(strcat('v',int2str([1:9]'))); %构造顶点字符串
G=graph(d,s); %构造无向图
T=minspantree(G); %用默认的Prim算法求最小生成树
L=sum(T.Edges.Weight) %计算最小生成树的权重 T.Edges %显示最小生成树的边集 h=plot(G,'Layout','circle','NodeFontSize',12); %画无向图
highlight(h,T,'EdgeColor','r','LineWidth',2) %加粗最小生成树的边
例7求最大流并标注
clc, clear, close all
a = zeros(6); a(1,[2,4])=[8,7];
a(2,[3,4])=[9,5]; a(3,[4,6])=[2,5];
a(4,5)=9; a(5,[3,6])=[6,10];
G = digraph(a);
H=plot(G,'EdgeLabel',G.Edges.Weight);
[M,F]=maxflow(G,1,6)
F.Edges, highlight(H,F, 'EdgeColor','r','LineWidth',1.5)
Matlab图论工具箱求解最大流的命令graphmaxflow,只能解 决权重都为正值,且两个顶点之间不能有两条弧的问题。若有两个顶点可添加一个虚拟顶点,例在顶点4和顶点3之间加入一个虚 拟顶点9,并添加两条弧,删除顶点4到顶点3的权重为2的弧,加 入的两条弧的容量都是2,不是1,因为这是容量,这一条需要2的容量,不是距离。
最小费用最大流问题就是要求一个从发点s到收点t的最大流,使流的总输送费用bf的和取最小值
function [flow,val]=mincostmaxflow(rongliang,cost,flowvalue);
%最小费用最大流函数
%第一个参数:容量矩阵;第二个参数:费用矩阵;
%前两个参数必须在不通路处置零
%第三个参数:指定容量值(可以不写, 表示求最小费用最大流)
%返回值flow 为可行流矩阵, val 为最小费用值
例8求最大流最小费用
global M num
c=zeros(6);u=zeros(6);
c(1,2)=2;c(1,4)=8;c(2,3)=2;c(2,4)=5;
c(3,4)=1;c(3,6)=6;c(4,5)=3;c(5,3)=4;c(5,6)=7;
u(1,2)=8;u(1,4)=7;u(2,3)=9;u(2,4)=5;
u(3,4)=2;u(3,6)=5;u(4,5)=9;u(5,3)=6;u(5,6)=10;
num=size(u,1);M=sum(sum(u))*num^2;
[f,val]=mincostmaxflow(u,c)
其中需要添加两个外置函数
其一
function path=floydpath(w)
global M num
w=w+((w==0)-eye(num))*M;
p=zeros(num);
for k=1:num
for i=1:num
for j=1:num
if w(i,j)>w(i,k)+w(k,j)
w(i,j)=w(i,k)+w(k,j);
p(i,j)=k;
end
end
end
end
if w(1,num)==M
path=[];
else
path=zeros(num);
s=1;
t=num; m=p(s,t);
while ~isempty(m)
if m(1)
s=[s,m(1)];t=[t,t(1)];t(1)=m(1);
m(1)=[];
m=[p(s(1),t(1)),m,p(s(end),t(end))];
else
path(s(1),t(1))=1;s(1)=[];m(1)=[];t(1)=[];
end
end
end
其二
function [flow,val]=mincostmaxflow(rongliang,cost,flowvalue);
global M
flow=zeros(size(rongliang));
allflow=sum(flow(1,:));
if nargin<3
flowvalue=M;
end
while allflow<flowvalue
w=(flow<rongliang).*cost-((flow>0).*cost)';
path=floydpath(w);
if isempty(path)
val=sum(sum(flow.*cost));
return;
end
theta=min(min(path.*(rongliang-flow)+(path.*(rongliang-flow)==0).*M));
theta=min([min(path'.*flow+(path'.*flow==0).*M),theta]);
flow=flow+(rongliang>0).*(path-path').*theta;
allflow=sum(flow(1,:));
end
val=sum(sum(flow.*cost));
两个新建函数,点击上面的运行(F9好像不行),自动命名为“函数名.m”,选择“添加到路径”,单独运行会提示输入参数不足,然后运行主文件
常用函数
旅行商问题:用图论的术语说,就是在一个赋权完全图(每对不同的顶点之间都恰连有一条边相连)中,找出一个有最小权的Hamilton圈。称这种圈为最优圈。目前较好的算法的是每次从已得到的圈中选两个点vi和vj,比较v(i)v(j)+v(i+1)v(j+1)与v(i)v(i+1)+v(j)v(j+1)进行迭代
例9旅行商
clc, clear, close all
global a
a=zeros(6);
a(1,2)=56;a(1,3)=35;a(1,4)=21;a(1,5)=51;a(1,6)=60;
a(2,3)=21;a(2,4)=57;a(2,5)=78;a(2,6)=70; a(3,4)=36;
a(3,5)=68;a(3,6)=68; a(4,5)=51;a(4,6)=61;a(5,6)=13;
a=a+a';L=size(a,1); %返回矩阵的行数,2则是列数
c1=[5 1:4 6 ]; %5/1:4/6
[circle,long]=modifycircle(c1,L); %
c2=[5 6 1:4 ]; %改变初始圈,该算法的最后一个顶点不动
[circle2,long2]=modifycircle(c2,L);
if long2<long
long=long2;
circle=circle2;
end
circle,long
外置函数
function [circle,long]=modifycircle(c1,L);
global a
flag=1;
while flag>0
flag=0;
for m=1:L-3
for n=m+2:L-1
if a(c1(m),c1(n))+a(c1(m+1),c1(n+1))<a(c1(m),c1(m+1))+a(c1(n),c1(n+1))
flag=1;
c1(m+1:n)=c1(n:-1:m+1);
end
end
end
end
long=a(c1(1),c1(L));
for i=1:L-1
long=long+a(c1(i),c1(i+1));
end
circle=c1;
若不使用图论的解法,还可以用线性0-1规划
例10
clc, clear
a=readmatrix('data4_15.xlsx');
a(isnan(a))=0; %把NaN替换为0
b=zeros(10); %邻接矩阵初始化
b([1:end-1],[2:end])=a; %邻接矩阵上三角元素赋值
b=b+b'; %构造完整的邻接矩阵
n=10;
b([1:n+1:end])=1000000; %对角线元素换为充分大
prob=optimproblem;
x=optimvar('x',n,n,'Type','integer','LowerBound',0,'UpperBound',1);
u=optimvar('u',n,'LowerBound',0) %序号变量
prob.Objective=sum(sum(b.*x));
prob.Constraints.con1=[sum(x,2)==1; sum(x,1)'==1; u(1)==0];con2 = [1<=u(2:end); u(2:end)<=14];
for i=1:n
for j=2:n
con2 = [con2; u(i)-u(j)+n*x(i,j)<=n-1];
end
end
prob.Constraints.con2 = con2;
[sol, fval, flag]=solve(prob)
xx=sol.x; [i,j]=find(xx);
fprintf('xij=1对应的行列位置如下:\n')
ij=[i'; j']
其中b([1:n+1:end])=1000000;的原理是矩阵b的元素按列排序记为b(n)如3*3的矩阵元素顺序就为