图论基础知识
图：$G=(V, E)$
点：$V={v_1,v_2,v_3...}$（结点数：n）
边：$E={e_1,e_2,e_3...}$（边数：m）
- 无向边-无向图$e_k$
- 有向边-有向图$e_k=(v_i,v_j)$ （$v_i$和$v_i$是相邻结点）
子图：给定图$G=(V, E)$，存在$G^=(V^, E^)$，满足$V^ \subseteq V$，$E^` \subseteq E$
支撑子图、生成子图：子图$G=(V, E)$中$V^`=V$
如何在计算机中表示一个图？
- 邻接矩阵：连通为1，否则为0 （可以表示子环，不可表示重边）
- 权矩阵：边的权值代替邻接矩阵的1
- 其他：关联矩阵、边列表、正向表、逆向表、邻接表
常见路的问题
道路、链：由一个结点到另一结点的路径
连通图：无向图中，任意两结点间都存在道路
判断两个结点间是否存在道路
- 广探法 (Breadth First Search,BFS)
- 深探法 (Depth First Search,DFS)
最短路径
- 单源最短路：两个结点间的最短路径 (Dijkstra)
- 某结点到其他各结点的最短路径 (Flord)
- 任意两结点间的最短路径 (Flord)
迪克斯特拉算法
%================================================================%
% 迪克斯特拉算法(贪心算法)：单源最短路径
% 输入：代价矩阵(M:不通)
% 输出：d(距离)，index1，index2
%================================================================%
clear all;
close all;
clc;
M = inf;
a(1,:)=[0,50,M,40,25,10];
a(2,:)=[zeros(1,2),15,20,M,25];
a(3,:)=[zeros(1,3),10,20,M];
a(4,:)=[zeros(1,4),10,25];
a(5,:)=[zeros(1,5),55];
a(6,:)=zeros(1,6);
% a(1,:)=[ 0, 50, M, 40, 25, 10];
% a(2,:)=[50, 0, 15, 20, M, 25];
% a(3,:)=[ M, 15, 0, 10, 20, M];
% a(4,:)=[40, 20, 10, 0, 10, 25];
% a(5,:)=[25, M, 20, 10, 0, 55];
% a(6,:)=[10, 25, M, 25, 55, 0];
%若输入的矩阵不是上三角矩阵，则语句 a=a+a'; 省略
a=a+a';
p_b(1:length(a)) = 0;
p_b(1) = 1;
index_1 = 1;
index_2 = ones(1,length(a));
d(1:length(a)) = inf;
d(1) = 0;
temp = 1;
while sum(p_b)<length(a)
tb = find(p_b==0);
d(tb) = min(d(tb),d(temp)+a(temp,tb));
tmpb = find(d(tb)==min(d(tb)));
temp = tb(tmpb(1));
p_b(temp) = 1;
index_1 = [index_1,temp];
index = index_1(find(d(index_1)==d(temp)-a(temp,index_1)));
if length(index)>=2
index = index(1);
end
index_2(temp) = index;
end
d, index_1, index_2
弗洛伊德算法
%================================================================%
% 弗洛伊德算法(动态规划)：两顶点之间的最短路径
% 输入：代价矩阵(M:不通)
% 输出：b(距离矩阵)，Path(路径矩阵)
%================================================================%
clear all;
close all;
clc;
M = inf;
a(1,:)=[0,50,M,40,25,10];
a(2,:)=[zeros(1,2),15,20,M,25];
a(3,:)=[zeros(1,3),10,20,M];
a(4,:)=[zeros(1,4),10,25];
a(5,:)=[zeros(1,5),55];
a(6,:)=zeros(1,6);
% a(1,:)=[ 0, 50, M, 40, 25, 10];
% a(2,:)=[50, 0, 15, 20, M, 25];
% a(3,:)=[ M, 15, 0, 10, 20, M];
% a(4,:)=[40, 20, 10, 0, 10, 25];
% a(5,:)=[25, M, 20, 10, 0, 55];
% a(6,:)=[10, 25, M, 25, 55, 0];
%若输入的矩阵不是上三角矩阵，则语句 b=a+a'; 省略
D=a+a';
Path = zeros(length(D));
for k=1:length(D)
for i=1:length(D)
for j=1:length(D)
if D(i,j)>D(i,k)+D(k,j)
D(i,j) = D(i,k)+D(k,j);
Path(i,j) = k;
end
end
end
end
D, Path
%求路径：floyd的后续指令
function pathway = road(path,v1,v2)
% 路径起点 辅助点 循环变量起点
pathway = v1;
q = v1;
k = 1;
while path(q,v2)~=v2 % q至v2后点不是v2
k = k+1;
pathway(k) = path(q,v2); % 路径增加点
q = path(q,v2); % 辅助点为新点,到终点v2结束
end
pathway(k+1) = v2;
旅行商问题(TSP)
哈密顿回路：无向连通图中，一条经过所有结点的初级(无重边、重点)回路。回路、圈：起始点相同的道路
概念： 重边：一对结点存在多条边（多重图）
kruskal.m（避圈法）
%================================================================%
% 最小生成树算法,通过kruskal算法求最优树,并给出相应图像
% 输入：S、E、W
% 输出：最小生成树的权
%================================================================%
clear all;
close all;
clc;
S=[1 1 2 2 3 4 4 5 5 6 6]; % 起始节点向量
E=[2 6 3 5 4 1 6 3 4 2 5]; % 终止节点向量
W = [.41 .29 .51 .32 .50 .45 .38 .32 .36 .29 .21]; % 边权值向量
DG = sparse(S,E,W);
% sparse：产生稀疏矩阵
% S：对应要形成矩阵的行位置
% E：对应要形成矩阵的列位置
% W：对应要形成矩阵对应位置的权值
% DG =
% (4,1) 0.4500
% (1,2) 0.4100
% (6,2) 0.2900
% (2,3) 0.5100
% (5,3) 0.3200
% (3,4) 0.5000
% (5,4) 0.3600
% (2,5) 0.3200
% (6,5) 0.2100
% (1,6) 0.2900
% (4,6) 0.3800
view(biograph(DG,[],'ShowArrows','off','ShowWeights','on'))
[ST,pred] = graphminspantree(DG) % 寻找最小生成树
view(biograph(ST,[],'ShowArrows','off','ShowWeights','on'))
W=sum(sum(ST)); % 最小生成树的权
disp(['最小生成树的权为：',num2str(W)]);
%================================================================%
% 最小生成树算法,通过kruskal算法求最优树
%================================================================%
clear all;
close all;
clc;
n = size(w,1);
% 求点边矩阵
k = 1;
for i=1:(n-1)
for j=(i+1):n
if w(i,j)~=inf
b(:,k)=[i;j;w(i,j)];
k=k+1;
end
end
end
m = size(b,2);
[b,I] = sortrows(b',3),b=b' % 排序
for i=1:m
if t(b(1,i))~=t(b(2,i)) % 不在同一子树
T(1:2,k) = b(1:2,i);
c = c+b(3,i);
tmin = min(t(b(1,i)),t(b(2,i)));
tmax = max(t(b(1,i)),t(b(2,i)));
for j=1:n % 更新子树的最小标号点
if t(j)==tmax
t(j) = tmin;
end
end
k = k+1;
end
if k==n
break;
end
end
MinTree_Prim.m
function [Tree,Weight,Weight1] = MinTree_Prim(M)
% 最小生成树的普利姆算法
% 输入：M 权值矩阵(一维)
% 输出：
% Tree ：最小生成树的边的集合
% Weight ：与选中边对应的权重
% Weight1 ：最小生成树的权之和
n = length(M); % 求图的顶点数
M(M==0) = inf;
k = 1:n;
listV(k) = 0;
listV(1) = 1;
e = 1;
while e<n
min = inf;
for i = 1:n
if listV(i) == 1
for j = 1:n
if listV(j) == 0 && min > M(i,j)
min = M(i,j);
B = M(i,j);
S = i;
D = j;
end
end
end
end
listV(D) = 1;
distance(e) = B;
source(e) = S;
destination(e) = D;
e = e+1;
end
Tree = [source;destination];
for g =1:e-1
Weight(g) = M(Tree(1,g),Tree(2,g));
end
Weight1 = sum(Weight);
Prim.m
clc;
clear;
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 length(result)~=length(a)-1
temp = a(p,tb);
temp = temp(:);
d = min(temp);
[jb,kb] = find(a(p,tb)==d);
j = p(jb(1));
k = tb(kb(1));
result = [result,[j;k;d]];
p = [p,k];
tb(find(tb==k)) = [];
end
result
% result 的第一、二、三行分别表示生成树边的起点、终点、权集合
MaxTree_Kruskal.m
function [Tree,Weight] = MaxTree_Kruskal(M, flag)
% 最大生成树的Kruskal算法，调用方法根据图的权值矩阵而定
% 当权值矩阵为三维时，
% 调用方法为：[Tree,Weight] = Ch19_Tree_Kruskal(M,1)，这里的1可以用任何数来代替
% 当“权值矩阵”为一维时，调用方法为：[Tree,Weight] = Ch19_Tree_Kruskal(M)
% 当“权值矩阵”为一维时：只需将求最小支撑树的权值矩阵中为inf的元素修改为-inf即可
% 输入
% M ：图（树）权值矩阵(3×n)，n为图的边数；每列的3个元素:边的起点、终点和权值。
% flag ：变量个数控制参数
% 输出
% Tree ：最大生成树的边的集合
% Weight ：最大生成树的权之和
% 注：图的边不重复编号
% 如果图的权值矩阵为1维，则进行如下转换
if nargin == 1
n = size(M,2); % 求图的顶点数
m = sum(sum(M~=0))/2; % 求图的边数
M1 = zeros(3,m);
k = 1;
for i = 1:n
for j = i+1:n
if M(i,j) ~= 0
M1(1,k) = i;
M1(2,k) = j;
M1(3,k) = M(i,j);
k = k+1;
end
end
end
else
M1 = M;
end
% 如果图的权值矩阵为3维，则直接进行求解
n = max(max(M1(1:2,:))); % 求图的顶点数
m = size(M1,2); % 求图的边数
B = -sortrows(-M1',3); % 按权的非增顺序重新排列边
B = B';
Tree = [];
Weight = 0;
k = 1;
q = 1:n;
for i = 1:m
if q(B(1,i)) ~= q(B(2,i))
Tree(1:2,k) = B(1:2,i);
Weight = Weight+B(3,i);
k = k+1;
qmin = min(q(B(1,i)),q(B(2,i)));
qmax = max(q(B(1,i)),q(B(2,i)));
for j = 1:n
if q(j) == qmax
q(j) = qmin;
end
end
end
if k == n
break;
end
end
MaxTree_Prim.m：
function [Tree,Weight,Weight1] = MaxTree_Prim(M)
% 最大生成树的Prim算法
% 输入:
% M ：图（树）权值矩阵(一维)：只将求最小支撑树的权值矩阵中的inf的元素修改为-inf即可
% 输出:
% Tree ：最大生成树的边的集合
% Weight ：与边对应的权值
% Weight1 ：最大生成树的权之和
n = length(M); % 求图的顶点数
M(M==0) = inf;
k = 1:n;
listV(k) = 0;
listV(1) = 1;
e = 1;
while e<n
max = 0; % 设定最大值的初值
for i = 1:n
if listV(i) == 1
for j = 1:n
if listV(j) == 0 && max < M(i,j)
max = M(i,j);
B = M(i,j);
S = i;
D = j;
end
end
end
end
listV(D) = 1;
distance(e) = B;
source(e) = S;
destination(e) = D;
e = e+1;
end
Tree = [source;destination];
for g =1:e-1
Weight(g) = M(Tree(1,g),Tree(2,g));
end
Weight1 = sum(Weight);
最佳匹配（二部图）
BestFit.m
%================================================================%
% 最佳匹配（二部图）
% 输入：m（行）
% n（列）
% 代价矩阵A（到达各顶点的权值）
% 输出：M矩阵
%================================================================%
clear all;
close all;
clc;
m = 5; % 上
n = 4; % 下
A = [ 4 5 5 1
2 2 4 6
4 2 3 3
5 0 2 1
6 3 1 4 ];
for i=1:n
L(i,1) = 0;
L(i,2) = 0;
end
% 初始化可行点标记 L
for i=1:n
for j=1:n
if L(i,1)<A(i,j)
L(i,1) = A(i,j);
end
R(i,j) = 0;
end
end
% 生成子图 G_l
for i=1:n
for j=1:n
if (L(i,1) + L(j,2))==A(i,j)
G_l(i,j) = 1;
else
G_l(i,j) = 0;
end
end
end
ii=0;
jj=0;
% 获得仅含 G_l 的一条边的初始匹配 M
for i=1:n
for j=1:n
if G_l(i,j)
ii = i;
jj = j;
break;
end
end
if ii
break;
end
end
R(ii,jj) = 1;
for i=1:n
S(i) = 0;
T(i) = 0;
Nis(i)=0;
end
while 1
for i=1:n
k = 1;
for j=1:n
if R(i,j)
k = 0;
break;
end
end
if k
break;
end
end
if k==0
break;
end % 获得最佳匹配 M, 算法终止
S(1) = i;
js_s = 1;
js_t = 0;
while 1
js_n = 0;
for i=1:js_s
for j=1:n
if G_l(S(i),j)
js_n = js_n+1;
Nis(js_n) = j;
for k=1:js_n-1
if Nis(k)==j
js_n = js_n-1;
end
end
end
end
end
if js_n == js_t
pd = 1; % 判断 NL(S)=T?
for j=1:js_n
if Nis(j)~=T(j)
pd = 0;
break;
end
end
end
if js_n==js_t & pd
al = Inf; % 如果 NL(S)=T, 计算 al, Inf 为∞
for i=1:js_s
for j=1:n
pd = 1;
for k=1:js_t
if T(k)==j
pd = 0;
break;
end
end
if pd & al>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:js_s
L(S(i),1) = L(S(i),1)-al;
end % 调整可行点标记
for j=1:js_t
L(T(j),2) = L(T(j),2)+al;
end % 调整可行点标记
for i=1:n
for j=1:n % 生成子图 G_l
if L(i,1)+L(j,2)==A(i,j)
G_l(i,j) = 1;
else
G_l(i,j) = 0;
end
R(i,j) = 0;
k = 0;
end
end
ii = 0;
jj = 0;
for i=1:n
for j=1:n
if G_l(i,j)
ii = i;
jj = j;
break;
end
end
if ii
break;
end
end % 获得仅含 G_l 的一条边的初始匹配 R
R(ii,jj) = 1;
break
else % NL(S)≠ T
for j=1:js_n
pd = 1; % 取 y∈ NL(S)\T
for k=1:js_t
if T(k)==Nis(j)
pd = 0;
break;
end
end
if pd
jj = j;
break;
end
end
pd=0; % 判断 y 是否为 M 的饱和点
for i=1:n
if R(i,Nis(jj))
pd = 1;
ii = i;
break;
end
end
if pd
js_s = js_s+1;
S(js_s) = ii;
js_t = js_t+1;
T(js_t) = Nis(jj); %S=S∪ {x}, T=T∪ {y}
else % 获得 G_l 的一条 M-增广路, 调整匹配 M
for k=1:js_t
R(S(k),T(k)) = 1;
R(S(k+1),T(k)) = 0;
end
if js_t==0
k = 0;
end
R(S(k+1),Nis(jj)) = 1;
break;
end
end
end
end
max_Weight = 0;
for i=1:n
for j=1:n
if R(i,j)
max_Weight = max_Weight+A(i,j);
end
end
end
R % 最佳匹配 M
max_Weight % 最佳匹配 M 的权, 程序结束
最大最小问题
Flord-Fulkerson最大流标号法
function minCost_maxFlow
clear;
clc;
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)
% 最小费用最大流函数
% 第一个参数：容量矩阵；第二个参数：费用矩阵；
% 前两个参数必须在不通路处置零
% 第三个参数：指定容量值（可以不写，表示求最小费用最大流）
% 返回值 flow 为可行流矩阵,val 为最小费用值
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));
Edmonds-Karp算法：
minCost_maxFlow.m
function [Tevent,Tactivity,KeyRoad,FinishTime] = NetTimeParamter(Matrix)
% 网络时间参数的计算,给出各时间和活动的时间参数，并给出关键线路
% 输入
% Matrix 网络图的赋权邻接矩阵
% 输出
% Tevent 事件的时间及关键事件
% Tactivity 活动的时间及关键活动、虚活动标示
% KeyRoad 关键线路，按照事件代号给出
% FinishTime 项目完工时间
m = size(Matrix); % 计算赋权网络图邻接矩阵的维度
% 计算各事件的时间参数
Tevent = zeros(m(1),5); % 存放时间的矩阵，包括关键事件及时差
Tevent(:,1) =1:m(1); % 第1列存放事件代号
% 计算事件的最早可能发生时间，逐列顺序计算
for j = 2:m(1) % 列标
M = zeros(j,1); % 存放前面所有节点到即将计算的节点的开始时间
for i = 1:j-1 % 行标
if Matrix(i,j)~=inf
M(i) = Matrix(i,j);
M(i) = Tevent(i,2)+M(i);
else
%M(i) = Tevent(i,2); % 这条语句不能有
end
end
Tevent(i+1,2) = max(M); % 第2列，存放事件的最早可能开始时间
end
% 计算事件的最迟必须发生时间，逐行倒序计算
Tevent(m(1),3) = Tevent(m(1),2);
for i = m(1)-1:-1:1 % 行标
M = zeros(m(1)-i,1);
for j = i+1:m(1) % 列标
if Matrix(i,j) ~= inf;
M(j-i) = Matrix(i,j);
M(j-i) = Tevent(j,3)-M(j-i);
else
M(j-i) = Tevent(m(1),3); % 保证逆向求差时，事件点能够取到正确值
end
end
Tevent(i,3) = min(M); % 第3列，存放事件的最迟必须发生时间
end
% 计算各时间的时差并查找关键路线
q = 1;
for p = 1:m(1)
Tevent(p,5) = Tevent(p,3)-Tevent(p,2); % 第4列，存放事件的时差
if Tevent(p,3)-Tevent(p,2) == 0
Tevent(p,4) = 1; % 关键事件
index(q) = p; % 关键线路上的事件代号
q = q+1; % 关键线路的指针
else
Tevent(p,4) = 0; % 非关键事件
end
end
KeyRoad = index(1:end);
% 计算项目的完工时间
FinishTime = Tevent(m(1),2);
% 计算各项活动的时间参数，逐列顺序计算
n = sum(sum(Matrix~=inf))-m(1); % 计算活动的数量=非无穷大元素个数减去零元素的个数
% 存放活动时间参数的矩阵，包括活动代号、4个时间参数、是否关键活动和总时差
Tactivity = zeros(n,9);
% 第1列，存放活动序号
% 第2列，存放活动的起始节点代号
% 第3列，存放活动的终止节点代号
% 第4列，存放活动的最早可能开始时间
% 第5列，存放活动的最迟必须开始时间
% 第6列，存放活动的最早可能完成时间
% 第7列，存放活动的最迟必须完成时间
% 第8列，表明活动是否为关键活动或虚活动，1表示是，0表示不是，2表示虚活动
% 第9列，存放活动的总时差
Tactivity(:,1) = 1:n;
h = 1;
for j = 2:m(1) % 列标
for i = 1:j-1 % 行标
if Matrix(i,j)~=inf;
Tactivity(h,2) = i;
Tactivity(h,3) = j;
Tactivity(h,4) = Tevent(i,2);
Tactivity(h,5) = Tactivity(h,4)+Matrix(i,j);
Tactivity(h,7) = Tevent(j,3);
Tactivity(h,6) = Tactivity(h,7)-Matrix(i,j);
if Tactivity(h,6)-Tactivity(h,7) == 0
Tactivity(h,8) = 2; % 虚活动
elseif Tactivity(h,6)-Tactivity(h,4) == 0
Tactivity(h,8) = 1; % 关键活动
else
Tactivity(h,8) = 0; % 非关键活动
end
Tactivity(h,9)=Tactivity(h,6)-Tactivity(h,4);
h = h+1; % 计数器累加1，控制存放位置
else
Tactivity(h,2:8) = 0;
end
end
end