图论模型和算法
Gary哥哥的哥哥的哥哥 2021.1.30
- 本博客聚焦于matlab解决与图论相关的数学建模问题
- 读者在阅读之前,最好对数据结构或者算法设计与分析的图论知识有所了解
图论基本概念
假设读者有相关基础,因此图论最基本的概念这里省略 相关知识可见此博客
说明:
- 接近中心度: v顶点到所有顶点的最短路径的总和的倒数,越大证明越“中心”
- 中间中心度(表示v的频率): s到t顶点的最短路径数量,分子为路径中经过了v顶点的最短路径数量
- 特征向量中心度:(求矩阵特征向量可得)
相关函数
Matlab中的相关函数
- 这里不对相关算法展开详细解释,相关问题可以在数据结构等课程中自行学习
graphallshortestpaths | 求图中所有顶点对之间的最短距离 |
---|---|
graphconnredcomp | 找无 (有) 向图的 (强/弱) 连通分支 |
graphisreddag | 测试有向图是否含有圈 |
graphisomorphism | 确定一个图是否有生成树 |
graphmaxflow | 计算有向图的最大流 |
graphminspantree | 在图中找最小生成树 |
graphpred2path | 把前驱顶点序列变成路径的顶点序列 |
graphshortestpath | 求指定一对顶点间的最短距离和路径 |
graphtopoorder | 执行有向无圈图的拓扑排序 |
graphtraverse | 求从一顶点出发, 所能遍历图中的顶点 |
- 问题抽象化:
- 在解决图论问题时候,应该先用矩阵表示出图
- 稀疏矩阵和满矩阵的转化函数:sparse & full
% a b c d e f
w = [0 0 0 0 0 0 % a
2 0 0 0 0 0 % b
3 6 0 0 0 0 % c
0 5 0 0 0 0 % d
0 3 1 1 0 0 % e
0 0 0 2 4 0]; % f
W = sparse(w);
[Dist,Path] = graphshortestpath(W, 1, 6,'Directed',0)
%output:
Dist =
7
Path =
1 3 5 4 6
网络分析工具箱
- 为了方便计算一些更加复杂的图论中的变量,我给出一个功能强大的工具箱
- 非Matlab官方的,是MIT开发的一个工具箱,地址(那个zip文件就是了)
函数名 | 功能 |
---|---|
degrees | 求图中所有顶点的度,入度和出度 |
ave_neighbor_deg | 求图中所有顶点的相邻顶点平均度 |
closeness | 求图中所有顶点的接近中心度 |
node_betweenness_faster | 求图中所有顶点的中间中心度 |
edge_betweenness | 求图中所有边的中间中心度 |
… | … |
案例
n = 10; %顶点总数
% 给每个人都编号
Andre = 1; Betty = 2; Carol = 3; Dave = 4; Ed = 5;
Fanny = 6; Garth = 7; Hale = 8; Ike = 9; Jane =10;
% 根据图构造邻接矩阵
A = zeros(10);
A(Andre, [Betty, Carol, Dave, Fanny]) = 1;% Andre跟哪几个人有边
A(Betty, [Andre, Dave, Ed, Garth]) = 1;
A(Carol, [Andre, Dave, Fanny]) = 1;
A( Dave, [Andre, Betty, Carol, Ed, Fanny, Garth]) = 1;
A( Ed, [Betty, Dave, Garth]) = 1;
A(Fanny, [Andre, Carol, Dave, Garth, Hale]) = 1;
A(Garth, [Betty, Dave, Ed, Fanny, Hale]) = 1;
A( Hale, [Fanny, Garth, Ike]) = 1;
A( Ike, [ Hale, Jane]) = 1;
A( Jane, [ Ike]) = 1;
Cd = degrees(A)' /(n-1) % 计算点度中心度并且标准化.
Cc = closeness(A)*(n-1) % 计算接近中心度并标准化
最小哈密尔顿路径问题
function [order, totdist] = minhamiltonpath(D)
%Step 1: Setup linear progrm
N = size(D, 1);
idxs = nchoosek(1:N, 2);
dist = D(sub2ind([N, N], idxs(:, 1), idxs(:, 2)));
lendist = length(dist);
%There need to be two "trips" atached to each "stop"
Aeq = spalloc(N, length(idxs), N*(N-1));
for i = 1:N
whichIdxs = idxs(:,1) == i | idxs(:,2) == i;
Aeq(i,whichIdxs) = 1;
end
beq = 2*ones(N, 1);
intcon = 1:lendist;
lb = zeros(lendist, 1);
ub = ones(lendist, 1);
opts = optimoptions('intlinprog','Display','off');
%Step 2: Solve linear program, adding inequality constraints until
%subtours are eliminated from found solution
x = intlinprog(dist,intcon,[],[],Aeq,beq,lb,ub,opts);
tours = detectSubtours(x, idxs);
ntours = length(tours);
A = spalloc(0, lendist, 0);
b = [];
while ntours > 1
fprintf(1, '%i subtours\n', ntours);
for i = 1:ntours
rowIdx = size(A, 1)+1;
touri = tours{i}; % Extract the current subtour
variations = nchoosek(1:length(touri), 2);
for j = 1:length(variations)
whichVar = (sum(idxs==touri(variations(j,1)),2)) & ...
(sum(idxs==touri(variations(j,2)),2));
A(rowIdx, whichVar) = 1;
end
b(rowIdx) = length(touri)-1;
end
x = intlinprog(dist,intcon,A,b,Aeq,beq,lb,ub,opts);
tours = detectSubtours(x, idxs);
ntours = length(tours);
end
totdist = dist'*x;
order = tours{:};
function subTours = detectSubtours(x,idxs)
% 返回子旅途元胞数组,即图中的子圈
x = round(x); % 纠正不确切的整数
r = find(x);
substuff = idxs(r,:); % 旅行线的节点对
unvisited = ones(length(r),1); % 跟踪未访问的旅行
curr = 1; % 正在评价的子旅途
startour = find(unvisited,1); % 第一个未访问的旅行
while ~isempty(startour)
home = substuff(startour,1);
nextpt = substuff(startour,2);
visited = nextpt;
unvisited(startour) = 0;
while nextpt ~= home
% 找以nextpt为起点的旅行
[srow,scol] = find(substuff == nextpt);
% 确定相应旅行的节点
trow = srow(srow ~= startour);
scol = 3-scol(trow == srow); % 1变2,2变1
startour = trow;
nextpt = substuff(startour,scol); % 子旅途的下一节点位置
visited = [visited,nextpt]; % 将节点加入子旅途
unvisited(startour) = 0; % 更新访问过的位置
end
subTours{curr} = visited; % 保存找到的子旅途
curr = curr + 1;
startour = find(unvisited,1);
end
end
tic
lat = [ 39.54 31.14 39.09 29.32 45.45 ...
43.52 41.50 40.49 38.02 37.52 ...
36.38 34.48 34.16 36.03 38.20 ...
36.38 43.48 31.51 32.02 30.14 ...
28.11 28.41 30.37 30.39 26.35 ...
26.05 23.08 20.02 22.48 25.00 ...
29.39 22.18 22.14 25.03 ];
lon = [116.28 121.29 117.11 106.32 126.41 ...
125.19 123.24 111.48 114.28 112.34 ...
117.00 113.42 108.54 103.49 106.16 ...
101.45 87.36 117.18 118.50 120.09 ...
113.00 115.52 114.21 104.05 106.42 ...
119.18 113.15 110.20 108.20 102.41 ...
90.08 114.10 113.35 121.31 ];
n = length(lat);
R = 6378.137;
dist = zeros(n);
for i = 1:n
for j = i+1:n
dist(i,j) = distance(lat(i),lon(i), lat(j),lon(j), R);
end
end
dist = dist + dist';
[order,totdist] = minhamiltonpath(dist);
plot(lon(order([1 : end ,1])), lat(order([1:end,1])),'o-')
灾情巡视路径问题 (CUMCM1998B)
-
三个县长把所有乡村跑一遍
-
兵分三路,总路程短一点,大家差不多时间一起停
-
分三个子图
-
分的方案很重要
-
先得出以o为顶点的最小生成树
-
再去分配
-
function [ncomp, icomp] = roadtree(ifplot) if nargin<1; ifplot = 'plot'; end [w, x,y, n, O] = road2graph('original', 'noplot'); [dist, path, pred] = graphshortestpath(w, O, 'directed', false); I = setdiff(1:n, O); J = pred(I); tree = zeros(n); for k = 1:length(I) i = I(k); j = J(k); tree(i,j) = w(i,j); tree(j,i) =w(j,i); end treeO = tree; treeO(:,O) = 0; treeO(O,:) = 0; %Find strongly or weakly connected components in graph [ncomp, icomp] = graphconncomp(sparse(treeO),'directed',false); if strcmp(ifplot,'plot') color = {'r','g','b','m','k','c'}; for k = 1:length(I) i = I(k); j = J(k); plot([x(i),x(j)], [y(i),y(j)]',['o-', color{icomp(i)}]) hold on end end
-
function [w, x,y, n, O] = road2graph(gtype, ifplot) if nargin<1; gtype = 'original'; ifplot = 'plot'; elseif nargin<2; ifplot = 'plot'; end nvil = 35; % number of villages ntown = 18; % number of towns n = nvil + ntown; ID = num2cell(nvil+1:nvil+ntown); % index of towns [A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R] = deal(ID{:}); %% coordinates of villages & towns for plot x = [61.7 51.3 52.0 51.3 44.4 35.9 33.3 35.5 22.9 20.9 ... % 1-10 18.9 9.4 9.1 1.6 0.0 8.2 13.9 15.1 21.1 26.1 ... % 11-20 24.7 20.3 28.6 29.9 30.5 43.5 44.1 50.6 59.0 60.8 ... % 21-30 64.4 68.9 68.1 75.1 78.5 67.2 65.2 58.4 45.9 28.6 ... % 31-35, A-E 19.9 13.1 1.9 8.1 14.9 20.9 26.5 39.3 34.9 56.3 ... % F- O 50.5 56.2 61.1 ]'; % P- R y = [28.9 23.6 19.0 6.9 24.6 25.4 19.2 13.3 12.7 0.0 ... % 1-10 19.2 11.8 24.3 26.7 39.7 47.2 45.9 35.6 28.3 33.6 ... % 11-20 39.8 48.6 46.7 54.2 37.4 41.7 47.7 45.5 45.7 53.7 ... % 21-30 45.8 50.5 40.5 36.5 42.7 34.2 25.9 20.9 15.9 15.3 ... % 31-35, A-E 8.8 17.6 17.6 37.7 28.6 41.2 28.7 32.6 44.3 29.6 ... % F- O 35.9 50.8 38.8 ]'; % P- R %% distance matrix w = zeros(n); w(1, A) = 10.3; w(1, B) = 5.9; w(1, C) = 11.2; w(1, O) = 6.0; w(2, 3) = 4.8; w(2, 5) = 8.3; w(2, O) = 9.2; w(3, C) = 7.9; w(3, D) = 8.2; w(4, 8) = 20.4; w(4, D) = 12.7; w(5, 6) = 9.7; w(5, D) = 11.3; w(5, M) = 11.4; w(6, 7) = 7.3; w(6, L) = 11.8; w(6, M) = 9.5; w(7, D) = 15.1; w(7, E) = 7.2; w(7, L) = 14.5; w(8, E) = 8.0; w(9, E) = 7.8; w(9, F) = 5.6; w(10,F) = 10.8; w(11,E) = 14.2; w(11,G) = 6.8; w(11,J) = 13.2; w(12,F) = 12.2; w(12,G) = 7.8; w(12,H) = 10.2; w(13,14)= 8.6; w(13,G) = 8.6; w(13,I) = 16.4; w(13,J) = 9.8; w(14,15)= 15.0; w(14,H) = 9.9; w(15,I) = 8.8; w(16,17)= 6.8; w(16,I) = 11.8; w(17,22)= 6.7; w(17,K) = 9.8; w(18,I) = 8.2; w(18,J) = 8.2; w(18,K) = 9.2; w(19,20)= 9.3; w(19,J) = 8.1; w(19,L) = 7.2; w(20,21)= 7.9; w(20,25)= 6.5; w(20,L) = 5.5; w(21,23)= 9.1; w(21,25)= 7.8; w(21,K) = 4.1; w(22,23)= 10.0; w(22,K) = 10.1; w(23,24)= 8.9; w(23,N) = 7.9; w(24,27)= 18.8; w(24,N) = 13.2; w(25,M) = 12.0; w(25,N) = 8.8; w(26,27)= 7.8; w(26,N) = 10.5; w(26,P) = 10.5; w(27,28)= 7.9; w(28,P) = 12.1; w(28,Q) = 8.3; w(29,P) = 15.2; w(29,Q) = 7.2; w(29,R) = 7.9; w(30,32)= 10.3; w(30,Q) = 7.7; w(31,32)= 8.1; w(31,33)= 7.3; w(31,R) = 9.2; w(32,33)= 19.0; w(32,35)= 14.9; w(33,35)= 20.3; w(33,A) = 7.4; w(34,35)= 8.2; w(34,A) = 11.5; w(34,B) = 17.6; w(A, B) = 12.2; w(A, R) = 8.8; w(B, C) = 11.0; w(I, J) = 15.8; w(C, O) = 11.5; w(M, N) = 14.2; w(M, O) = 19.8; w(O, P) = 10.1; w(O, R) = 12.9; w = sparse(w'); % upper triangular matrix to lower trianglular matrix %% plot graph if strcmp(ifplot,'plot') [i,j] = find(w>0); plot([x(i),x(j)]', [y(i),y(j)]','o:b','linewidth',4); end %% if strcmp(gtype,'complete') w = graphallshortestpaths(w, 'directed', false); end
-
-
总路程:三个子图最小哈密尔顿回路最小
-
均衡化 最小哈密尔顿回路三个中max-min小
-
具体过程见上面的pdf和代码
-