数模培训第二周——图论模型

图论中最短路算法与程序实现

  图论中的最短路问题(包括无向图和有向图)是一个基本且常见的问题。主要的算法有Dijkstra算法和Floyd算法。

Floyd算法

简介

  Floyd‐Warshall算法（英语：Floyd‐Warshall algorithm或简写为Floyd algorithm），中文亦称弗洛伊德算法，是解决任意两点间的最短路径的一种算法，可以正确处理无向图或有向图（可以有负权重，但不可存在负权回路）的最短路径问题。Floyd算法与迪杰斯特拉算法或贝尔曼福特算法相比，能够一次性的求出任意两点之间的最短路径，后两种算法运行一次只能计算出给定的起点和终点之间的最短路径。当然，Floyd算法计算的时间也要高于后两种算法，其算法核心的步骤由三层循环构成。

算法复杂度

时间复杂度$O(n^3)$，$n$为点的个数

引入

  一些城市之间有公路，有些则没有，如下图。为了方便安排行程以及节省路费，游客需要提前知道任意两个城市之间的最短路径和路程。

  图中有4个城市，每天边上的数字代表了该段公路的路程，注意这里的公路都是单向的。在求解最短路径之前，我们需要用邻接矩阵D来存储这个图。

  比如$A1$到$A2$的路程为$35$，则$D(1,2)$为$35$；$A2$到$A1$的路程为$36$，则$D(2,1)$为$36$；$A1$无法到达$A4$，则设置$D(1,4)$为无穷大$\infty$；$A1$到他自己的距离为$0$，所以$D(1,1)$为$0$。接着用路由矩阵$R$来存储两点之间的最短路径。
  首先，如果任意两点之间不能经过第三个点时，很明显这些城市之间的最短路程就是初始路程，两个矩阵都不变。
  接着，假设任意两点之间只可以经过
$A1$，该怎么求任意两点之间的最短路径呢？只需判断$D(i,1)+D(1,j)$是否比$D(i,j)$小即可。$D(i,j)$表示从顶点i到顶点j的路程，$D(i,1)+D(1,j)$表示从i号顶点到$1$号顶点，再从$1$号顶点到$j$号顶点的路程之和。其中$i$和$j$都是$1\sim n$循环，$n$为顶点数，即图上的城市数量之和。

  例如原本$A4$无法直接到达$A2$，故$D(4,2)$为无穷大；现在可以从$A4$经过$A1$到达$A2$，路程之和为$105$，即符合条件：$\infty =D\left(4,2\right)>D\left(4,1\right)+D\left(1,2\right)=105$，因此我们可以对两个矩阵进行更新： $D\left(4,2\right)=105,R\left(4,2\right)=R\left(4,1\right)=1$， 。这里$R(4,2)=1$表示从顶点$4$到顶点$2$的最短路径是先经过顶点$1$。

  同理， $\infty =D\left(4,3\right)>D\left(4,1\right)+D\left(1,3\right)=155$，即原本$A4\to A3=\infty$，现在经过$A1$中转后 ，于是$A4$到$A3$的最短路程变为$155$，$D(4,3)=155，R(4,3)=R(4,1)=1$

现在假设在允许经过$A1$中转的基础上，任意两点之间还可以经过$A2$中转，于是有：

原本：$A1\to A3=85$，现在：$A1\to A2\to A3=75$，满足$85=D\left(1,3\right)>D\left(1,2\right)+D\left(2,3\right)=75$，故更新矩阵：$D(1,3)=75,R(1,3)=R(1,2)=2$

原本：$A3\to A1=90$，现在：$A3\to A2\to A1=77$，满足$90=D\left(3,1\right)>D\left(3,2\right)+D\left(2,1\right)=77$$，故更新矩阵：$D(3,1) = 77, R(3,1) = R(3,2) = 2$

原本：$A4\to A1\to A3=155$，现在：$A4\to A1\to A2\to A3=145$，满足$155=D\left(4,3\right)>D\left(4,1\right)+D\left(1,3\right)=145$，故更新矩阵：$D(4,3)=145,R(4,3)=R(4,1)=1$

同理，在允许经过$A1$和$A2$中转的基础上，任意两点之间还可以经过$A3$中转，按照上文的方法，可以将邻接矩阵$D$和路由矩阵$R$更新为：

最后，允许经过全部顶点中转，可以得到最终的邻接矩阵$D$和路由矩阵$R$：

此时邻接矩阵$D$就是任意两点之间的最短路程矩阵，路由矩阵存储了对应的最短路径。举个栗子：求$A4$到$A3$的最短路径。

最短路程：$D(4,3)=145$；
最短路径：$R(4,3)=1,R(1,3)=2,R(2,3)=3$，故最短路径为$A4\to A1\to A2\to A3$

例题

已知50个点之间相互连接信息见表求最短距离矩阵，最短路径。

50个点之间相互连接信息矩阵

Floyd.m

function [distance, route] = Floyd(A)
    % Floyd算法求最短路径
    % 输入：邻接矩阵x，数值型
    % 输出：最短路径的距离矩阵distance，路由矩阵route

    n = size(A, 1);
    distance = A;

    % 路由矩阵route记录最短路径
    for i = 1 : n
       for j = 1 : n
          route(i, j) = j; 
       end
    end

    % Floyd算法
    for k = 1 : n
        for i = 1 : n
            for j = 1 : n
                if distance(i, j) > distance(i, k) + distance(k, j)
                    distance(i, j) = distance(i, k) + distance(k, j);
                    route(i, j) = route(i, k);
                end
            end
        end
    end
end

displayPath.m

function path = displayPath(route, start, dest)
    % 打印出任意两点之间的最短路径
    % route : 路由矩阵
    % start : 起点index
    % dest : 终点index
    % path：最短路径
    path = [start];
    while 1
       if route(start, dest) ~= dest 
          start = route(start, dest);
          path = [path, start];
       else
           path = [path, dest];
           break;
       end
    end
end

main.m

%%
clear, close, clc

%% 初始距离矩阵
n = 50;
A = (ones(n, n) - eye(n, n)) * 100000;
A(1,2)=400;A(1,3)=450; A(2,4)=300;A(2,21)=230; A(2,47)=140;A(3,4)=600;
A(4,5)=210;A(4,19)=310;A(5,6)=230;A(5,7)=200; A(6,7)=320; A(6,8)=340;
A(7,8)=170;A(7,18)=160;A(8,9)=200;A(8,15)=285; A(9,10)=180; A(10,11)=150;
A(10,15)=160; A(11,12)=140; A(11,14)=130; A(12,13)=200; A(13,34)=400;
A(14,15)=190;A(14,26)=190; A(15,16)=170; A(15,17)=250; A(16,17)=140;
A(16,18)=130; A(17,27)=240; A(18,19)=204; A(18,25)=180; A(19,20)=140;
A(19,24)=175; A(20,21)=180; A(20,24)=190; A(21,22)=300; A(21,23)=270;
A(21,47)=350;A(22,44)=160;A(22,45)=270;A(22,48)=180;A(23,24)=240;
A(23,29)=210;A(23,30)=290;A(23,44)=150;A(24,25)=170;A(24,28)=130;
A(26,27)=140;A(26,34)=320;A(27,28)=190;A(28,29)=260;A(29,31)=190;
A(30,31)=240;A(30,42)=130;A(30,43)=210;A(31,32)=230;A(31,36)=260;
A(31,50)=210;A(32,33)=190;A(32,35)=140;A(32,36)=240;A(33,34)=210;
A(35,37)=160;A(36,39)=180;A(36,40)=190;A(37,38)=135;A(38,39)=130;
A(39,41)=310;A(40,41)=140;A(40,50)=190;A(42,50)=200;A(43,44)=260;
A(43,45)=210;A(45,46)=240;A(46,48)=280;A(48,49)=200;
% 矩阵对称
for i = 1 : n
    for j = 1 : i - 1
        A(i, j) = A(j, i);
    end
end

%% 输出距离矩阵到文件
fid = fopen('primaryDistance.txt', 'w');
for i = 1 : n
    for j = 1 : n
        fprintf(fid, '%4d ', A(i,j));
    end
    fprintf(fid, "\n");
end
fclose(fid);

%%
clear, close, clc

%% 导入初始距离矩阵
load primaryDistance.txt
A = primaryDistance;
clear primaryDistance
n = size(A, 1);

%% 调用Floyd.m算法计算最短距离矩阵以及路由矩阵
[distance, route] = Floyd(A);

%% 将最短距离矩阵写入文件
fid=fopen('shortDistance.txt', 'w');
for i = 1 : n
    for j =1 : n
        fprintf(fid, '%4d ', distance(i,j));
    end
    fprintf(fid, "\n");
end
fclose(fid);

%% 将最短路径写入文件
fid=fopen('shortPath.txt', 'w');
for i = 1 : n
    for j = 1 : i
        if i == j
            continue;
        else
            path = displayPath(route, i, j);
            pathNumber = length(path);
            str = sprintf('(%d 到 %d 的最短路径长度为：%d) %d', i, j, distance(i, j), path(1));
            fprintf(fid, '%s', str);
            for k = 2 : pathNumber
                strTmp = sprintf(' -> %d', path(k));
                fprintf(fid, '%s', strTmp);
            end
            fprintf(fid, '\n');
        end
    end
end
fclose(fid);

部分结果

shortPath
shortDistance

---

TSP问题

简介

  巡回旅行商问题(Traveling Salesman Problem，TSP)，也称为货郎担问题。最早可以追溯到1759年Euler提出的骑士旅行问题。该问题可简单描述为走遍n个城市的最短路。1948年，由美国兰德公司推动，TSP成为近代组合优化领域的一个典型难题。它已经被证明属于NP难题。该问题在数学建模竞赛中也经常碰到。
这类问题我们可以用计算机模拟，智能算法解决，这里我们介绍利用LINGO软件进行求解的方法。

问题1

  设有一个售货员从10个城市中的某一个城市出发，去其它9个城市推销产品。10个城市相互距离如下表。要求每个城市到达一次仅一次后，回到原出发城市。问他应如何选择旅行路线，使总路程最短。

方案比较

设有6个城市，下面矩阵代表了一种方案:

方案1

合格方案

1->2->3->4->5->6->1

$$\left[\begin{array}{llllll}0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ 1& 0& 0& 0& 0& 0\end{array}\right]$$

方案2

不合格方案

1->2->3->1 4->5->6->4

$$\left[\begin{array}{llllll}0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 1& 0& 0\end{array}\right]$$

模型建立

设城市之间距离用矩阵$d$来表示，$d_{ij}$表示城市$i$与城市$j$之间的距离
设0-1矩阵$X$用来表示经过各个城市之间的路线。设:

$$x_{ij}=\{\begin{array}{cc}0& \text{ 若城市 }i\text{ 不到城市 }j\\ 1& \text{ 若城市 }i\text{ 到城市 }j\end{array}$$

考虑每一个城市后只有一个城市，则有：
$$\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^n{x}_{ij}=1,i=1,\cdots ,n$$

考虑每个城市前只有一个城市，则有：

$$\sum _{\genfrac{}{}{0ex}{}{i=1}{i\ne j}}^n{x}_{ij}=1,i=1,\cdots ,n$$

---

但是该约束不能产生子圈。为此引入额外变量$u_i$，附加以下约束：

$$u_i-u_j+n{x}_{ij}\le n-1,1<i\ne j\le n$$

解释：
如$i$与$j$不会构成回路，若构成回路，有：
$$x_{ij}=1，x_{ji}=1$$

则：
$$u_i-u_j⩽-1，u_j-u_i⩽-1，从而0⩽-2，导致矛盾。$$

如$i，j$与$k$不会构成回路，若构成回路，则有：

$$x_{ij}=1，x_{jk}=1，x_{ki}=1$$

则：
$$u_i-u_j\le -1,u_j-u_k\le -1,u_k-u_i\le -1\text{ 从而 }0\le -3$$

其他情况一次类推。

总线性规划模型为：

$$\begin{gathered} \min z=\sum _{i=1}^n\sum _{j=1}^n{d}_{ij}{x}_{ij} \\ s.t\{\begin{array}{l}{\sum }_{\begin{array}{c}i=1\\ i\ne j\end{array}}^n{x}_{ij}=1,j=1,\cdots ,n\\ {\sum }_{\begin{array}{c}j=1\\ j\ne i\end{array}}^n{x}_{ij}=1,j=1,\cdots ,n\\ u_i-u_j+n{x}_{ij}\le n-1,1<i\ne j\le n\\ x_{\begin{array}{c}ij\end{array}}=0or1,\left(i,j=1,\cdots ,n\right)\\ u_i\text{为实数},\left(i=1,\cdots ,n\right)\end{array} \end{gathered}$$

模型求解

Model :

Sets :
row/1..10/ : u;
link(row, row) : d, x;

Endsets

Data :
d = 0 7 4 5 8 6 12 13 11 18,
    7 0 3 10 9 14 5 14 17 17,
    4 3 0 5 9 10 21 8 27 12,
    5 10 5 0 14 9 10 9 23 16,
    8 9 9 14 0 7 8 7 20 19,
    6 14 10 9 7 0 13 5 25 13,
    12 5 21 10 8 13 0 23 21 18,
    13 14 8 9 7 5 23 0 18 12,
    11 17 27 23 20 25 21 18 0 16,
    18 17 12 16 19 13 18 12 16 0;
@text('2.txt') = @writefor(link(i, j) | x(i, j) #gt# 0 : @name(x), ' ', x(i, j), @newline(1));

Enddata

Min = @sum(link : d * x);

@for(row(j) : @sum(row(i) | i #ne# j : x(i, j)) = 1);
@for(row(i) : @sum(row(j) | j #ne# i : x(i, j)) = 1);
@for(link(i, j) | i #ne# j #and# i #gt# 1 : u(i) - u(j) + 10 * x(i, j) <= 9);
@for(link : @bin(x));
@for(row : @free(u));

End
  Global optimal solution found.
  Objective value:                              77.00000
  Objective bound:                              77.00000
  Infeasibilities:                             0.1776357E-14
  Extended solver steps:                               0
  Total solver iterations:                           486
X( 1, 9)
X( 2, 3)
X( 3, 4)
X( 4, 1)
X( 5, 7)
X( 6, 5)
X( 7, 2)
X( 8, 6)
X( 9, 10)
X( 10, 8)
1->9->10->8->6->5->7->2->3->4->1

---

最优树问题

简介

树：连通且不含圈的无向图称为树．常用$T$表示。树中的边称为树枝，树中度为1的顶点称为树叶。

生成树：若$T$是包含图$G$的全部顶点的子图,它又是树,则称$T$是$G$的生成树.

最小生成树：设$T=(V,E_1)$是赋权图$G=(V,E)$的一棵生成树，称$T$中全部边上的权数之和为生成树的权，记为$w(T)$，即$w(T)=\sum w(e)$。如果生成树$T^{\ast }$的权$w(T^{\ast })$是$G$的所有生成树的权中最小值，则成$T^{\ast }$是$G$的最优树，即$w(T^{\ast })=\sum \min \left\{w(T)\right\}$

  在许多实际问题中，如在许多城市间建立公路网、输电网或通信网络，都可以归结为赋权图的最优树问题。如在一个城市中，对若干个居民点要供应自来水，已经预算出连接各点间管道的造价，要求给出一个总造价最小的铺设方案。

图论中最优树的的求解通常有两种算法：$Kruskal$算法（或避圈法）和$Prim$算法（破圈法）

例题

  某有10个城镇见下图，它们之间的距离见表。城镇1处有一条河流，现需要从各城镇之间铺设管道，使城镇1处的水可以输送到各城镇，求铺设管道最少的设计方式。

模型建立

设无向图共有$n$个节点，其赋权图的邻接矩阵为$d_{n\times n}$。$d_{ij}$表示节点$i$到$j$的距离。$d$为对称矩阵。令$d_{ii}=0$。现求根节点1到各节点生成的最优树，要求各线路上的权值和最小，其线性规划模型为：
决策变量：设：
$$x_{ij}=\{\begin{array}{ll}1& \text{ 节点 }i\text{ 与节点 }j\text{ 连通 }\\ 0& \text{ 节点 }i\text{ 与节点 }j\text{ 不连通 }\end{array}$$

目标函数为寻找一条从起始点1到各节点生成的最优树，要求各线路上的权值和最小，故目标函数为:

$$\min Z=\sum _{i=1}^n\sum _{j=1}^n{d}_{ij}\cdot x_{ij}$$

• 对起始点1知道有一条路出去，故有：${\sum }_{j=2}^n{x}_{1j}⩾1$
• 对其余各节点，恰有一条路进入，有：
  $$\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^n{x}_{ji}=1i=2,3,\cdots ,n$$
• 所有节点不出现圈，约束为：$u_i-u_j+n{x}_{ij}⩽n-1i,j=1,2,\cdots ,n$

故总线性规划模型为：

$$\begin{gathered} \min Z=\sum _{i=1}^n\sum _{j=1}^n{d}_{ij}\cdot x_{ij} \\ s.t\{\begin{array}{l}{\sum }_{j=2}^n{x}_{1j}⩾1\\ {\sum }_{\begin{array}{c}j=1\\ j\ne i\end{array}}^n{x}_{ji}=1i=2,3,\cdots ,n\\ u_i-u_j+n{x}_{ij}⩽n-1i,j=1,2,\cdots ,n\\ x_{ij}=0or1\end{array} \end{gathered}$$

模型求解

Model :

Sets :
row/1..10/ : u;
link(row, row) : d, x;

Endsets

Data :
d = 0 8 5 9 12 14 12 16 17 22
    8 0 9 15 16 8 11 18 14 22
    5 9 0 7 9  11 7 12 12 17
    9 15 7 0 3 17 10 7 15 15
    12 16 9 3 0 8 10 6 15 15
    14 8 11 17 8 0 9 14 8 16
    12 11 7 10 10 9 0 8 6 11
    16 18 12 7 6 14 8 0 11 11
    17 14 12 25 15 8 6 11 0 10
    22 22 17 15 15 16 11 11 10 0;

@text('2.txt') = @writefor(link(i, j) | x(i, j) #gt# 0 : @name(x), '=', x(i, j), @newline(1));

Enddata

Min = @sum(link(i, j) | i #ne# j : d(i, j) * x(i, j));

@sum(row(j) | j #gt# 1 : x(1, j)) >= 1;
@for(row(i) | i #ne# 1 : @sum(row(j) | j #ne# i : x(j, i)) = 1);
@for(link : @bin(x));
@for(link(i, j) | i #ne# j : u(i) - u(j) + 10 * x(i, j) <= 9);

End
  Global optimal solution found.
  Objective value:                              60.00000
  Objective bound:                              60.00000
  Infeasibilities:                              0.000000
  Extended solver steps:                               4
  Total solver iterations:                           163
  X( 1, 2)=1
X( 1, 3)=1
X( 3, 4)=1
X( 3, 7)=1
X( 4, 5)=1
X( 5, 6)=1
X( 5, 8)=1
X( 7, 9)=1
X( 9, 10)=1

故最优树(最佳铺设管道方式)如下图：

---

如有错误以及可以改进的地方欢迎在下方评论区留言！