浙江理工大学2019数学建模校赛B题记录

Posted on | | Visitors:
Words count in article: 4,702 字 | Reading time ≈ 18 分

第一题

数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
w=[...
    0,   80, 150,   90, 140, 100, 120; 
120,     0, 180, 100,   60,   70, 110; 
100, 160,     0,   80,   90,   50,   70;
  inf,   inf,   inf,     0,   40,   70, 180; 
  inf,   inf,   inf,   30,     0,   80, 150; 
  inf,   inf,   inf,   80,   90,     0,   50;
  inf,   inf,   inf, 130, 140,   60,     0; ];
% 设甲、乙、丙三个矿区为顶点v1，v2，v3，
% A、B、C、D 三个冶炼厂为v4、v5、v6、v7。
% 该矩阵表示7个顶点间的距离
% 由于不是对称矩阵，第一行第二行的元素表示甲到乙还是乙到甲是有区别的。题目不严谨，而老师PPT里是前者。有待区分(TODO)
% TODO:graphallshortestpaths函数中的图参数 inf用0表示也可以？
% 该矩阵与题给数据是否完全一样，有待比较(TODO)

建立模型

构造有向图，G=(V,E, W)，顶点、边、权重。甲、乙、丙三个矿区为顶点v1，v2，v3，
% A、B、C、D 三个冶炼厂为v4、v5、v6、v7。TODO:这部分可以参考骆桦PPT或者相关书籍与论文。

xij>0，肯定大于0。即负数用表示同样含义的正数表示（TODO：怎么说呢）

由于可以转运：

IMPORTANT：为什么不是用该矩阵中的值作为顶点之间的距离，而是根据此矩阵计算最短路径？（因为题中说可以转运，而最短路径即包括了转运的含义）

问题中最优的定义：

矿区的矿物产量刚好达到（即等于）冶炼厂矿物需求量（TODO：矿物产量和矿物需求量这两个词是否合适，应参照题目或统一规定用词），这是两个约束条件。

IMPORTANT：因为题目是求矿区到各冶炼厂间矿石的最优调运方案，所以从最短路径中提取出两个矿区（甲乙）到四个冶炼厂（ABCD）的最短距离

分别用i =1, 2表示甲乙两个矿区，j =1，2，3，4表示A、B、C、D三个电厂， cij表示第i个矿区到第j个冶炼厂的最短距离，xij表示第i个矿区到第j个冶炼厂的调运量，ai表示第i个矿区的产量，bj表示第j个冶炼厂的需求量。

这里是产量和需求平衡的运输问题

一些假设

  • 花费与路径长度和运送量相关，所以假设花费=路径长度*运送量。在这一点可以尝试讨论，优化？
  • 假设调运方案中从矿区到各冶炼厂间调运吨数为整数，试试改成小数？

matlab程序

计算出的最短路径矩阵，前三行后四列是我们的结果，如下。

1
2
3
90   130   100   120
90     60     70   110
80     90     50     70
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
clc;clear;

w=[...
    0,   80, 150,   90, 140, 100, 120; 
120,     0, 180, 100,   60,   70, 110; 
100, 160,     0,   80,   90,   50,   70;
  0,   0,   0,     0,   40,   70, 180; 
  0,   0,   0,   30,     0,   80, 150; 
  0,   0,   0,   80,   90,     0,   50;
  0,   0,   0, 130, 140,   60,     0; ];

W = sparse(w);
d = graphallshortestpaths(W)
NodeIDs = {'甲', '乙','丙', 'A', 'B', 'C', 'D'}; % 结点标签，也就是h.Nodes(i).ID属性值
h = view(biograph(w, NodeIDs, 'ShowWeights', 'on'))
set(h.Nodes, 'shape', 'circle'); % 顶点画成圆形
h.EdgeType = 'segmented'; % 边的连接为线段
h.LayoutType = 'radial';
try
dolayout(h) % 刷新图形 会报错，故放在try块中，不影响结果，待深究
catch exception
end
h2 = view(biograph(d, NodeIDs, 'ShowWeights', 'on'));
h2.EdgeType =  'segmented'; % 边的连接为线段
h2.LayoutType = 'equilibrium';
try
dolayout(h2) % 刷新图形 会报错，故放在try块中，不影响结果，待深究
catch exception
end

TODO:matlab的图不怎么样。修改matlab参数、图论工具箱、网络分析工具箱、NetworkX

lingo程序

矩阵x和Objective Value是我们的最终结果：各矿区到各冶炼厂最优调度方案。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
model:
sets:
kuang/1..3 /:a; !甲乙丙矿区的产量;
chang/1..4 /:b; !ABCD四个冶炼厂的矿石用量;
link(kuang,chang):c,x;
endsets
data:
a = 700 500 500;
b = 400 300 400 600;
c = 90   130   100   120
	 90     60     70   110
	 80     90     50     70;
enddata
min = @sum(link: c * x);
@for(kuang(i):@sum(chang(j):x(i, j)) = a(i));
@for(chang(j):@sum(kuang(i):x(i, j)) = b(j));
end

TODO：检验答案正确性

135000.0

1
2
3
4
结果：三个矿区到四个冶炼厂的调度方案
400 0  200 100
0 300 200 0
0 0 0 500

待解决主要问题

图画的不行、论文（建立模型部分、假设部分），可能这只是第一题，不需要花太大精力。

第二题

几乎不用假设和建模，而且老师也给了两个lingo程序。但现在问题是老师给的第二个lingo程序看不太懂，不知道怎么把B题里的数据套进去。

我借的matlab书上都有直接给代码，但很长。

最大流

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
model:
sets:
nodes/s,1,2,3,t/;
arcs(nodes,nodes):c,f;
endsets
data:
c=0;
@text('fdata.txt')=f;
enddata
calc:
c(1,2)=20;c(1,3)=16;
c(2,4)=4;c(2,5)=14;
c(3,2)=10;c(3,4)=20;
c(4,5)=8;
endcalc
n=@size(nodes);
max=flow;
@for(nodes(i)|i #ne# 1 #and# i #ne#n:@sum(nodes(j):f(i,j))=@sum(nodes(j):f(j,i)));    
@sum(nodes(i):f(1,i))=flow;
@sum(nodes(i):f(i,n))=flow;
@for(arcs:@bnd(0,f,c));
end

最小费用最大流

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
model:
sets:
nodes/s,1,2,3,t/:d;
arcs(nodes,nodes):b,c,f;
endsets
data:
d=22 0 0 0 -22;
b=0;c=0;
enddata
calc:
b(1,2)=12;b(1,3)=3;
b(2,4)=18;b(2,5)=3;
b(3,2)=6;b(3,4)=9;
b(4,5)=6;
c(1,2)=20;c(1,3)=16;
c(2,4)=4;c(2,5)=14;
c(3,2)=10;c(3,4)=20;
c(4,5)=8;
endcalc
min=@sum(arcs:b*f);
@for(nodes(i):@sum(nodes(j):f(i,j))-@sum(nodes(j):f(j,i))=d(i));   
@for(arcs:@bnd(0,f,c));
End

第三题

可选方法

数学实验13-树算法08PPT里的后三种方法：

改进的贪心

模拟退火

修改过的prim

英文论文阅读

s commentary the outstanding Steiner tree papers

准确解

  • 暴力法

直接暴力，890+个可能位置。

Hanan的一个理论，被三个队伍引用，说只需要考虑和一个已知站x坐标相同并且和另外一个已知站y坐标相同的steiner点。可能的steiner点位置减少到63个。

在此基础上可以，删除角落里的点，可能的位置减少到31个。

  • 动态规划

虽然它还是要计算已知点的所有子集，但它避免测试所有steiner点的子集。

通过该这样的方法，几个小时可以找到50个可能位置，但似乎没有人这样做。

电脑芯片设计者通常喜欢近似解，因为网络的长度并不总是操作速度的限制条件。

近似解

一个近似算法是不用虚设点，直接使用最小生成树，Hwang证明这个长度不超过steiner树的150%，但这当然不满足题目的要求。

某队找到一个最好的steiner点，然后把它加进已知站，之后重复这个策略直到加了7个点或者不可能优化了。这个方法找到了一个可能解，但并不是每次都可以这样。

某队用了多种不同的启发式算法，并且用另外四个样例对他们的方法进行了测试。另外，他们在评估他们的算法时做了一个好的决定。关于最小生成树他们还给了一个优化方案，因为他们通常不能计算准确的steiner树。他们的启发式算法都用了贪心策略。他们想出了一种改进的克鲁斯卡尔算法，模仿克鲁斯卡尔但如果有利的话会使用steiner点。这个启发式算法找到了一种可能解。虽然这个算法在速度上优于上一队，但面对更大规模的问题时它不一定能找到相当好的解。

某队(也就是下边的这个论文)用了模拟退火，模拟退火根据当地重新排列规则随机从某种解移向某种解。移向一个解的概率取决于这两个解的花费和一个叫做温度的控制变量。通过一个合适的冷却速率，这个启发式算法最终停在一个几乎最优解。他们模拟了100次这个程序都得到了最优解。

总结

最后简单总结一下，有两种方法：

  • exact solution

    求确切的最优解，有两种方法。

    • 暴力枚举法

      通过减少steiner点可能的位置来节约时间，仍然是在用这种方法。

    • 动态规划

      似乎没有人用这个。

  • 近似解

    通过启发式算法，求近似最优解，有几种方法。

    • 单点贪心

      每次先找到一个最好的steiner点，将其加入已知站，直至N-2或不可能有更好的解了。

    • 改进克鲁斯卡尔算法

      模仿克鲁斯卡尔，但考虑了steiner点。

    • 模拟退火，这是最好的一个方法

      还未具体去看那篇论文(1991 B Finding optimal Steiner trees)

1991 B Finding optimal Steiner trees

模拟退火

过程

这是一种更有效的寻找所有可行的虚设点和已知点组合的方式。

我们从一个已给的虚设点和已知点组成的树开始，并且允许模拟退火程序创造一个新的配置。对每一个新配置来说，程序决定了最小生成树并且计算其长度。然后模拟退火程序在冷却计划的基础上决定是使用还是拒绝新配置。修改一个配置和评价这个修改方案都很简单。

这个方法很有普适性：任何初始配置都可能被使用，并且算法并不直接取决于固定的站点设置。

和暴力相反，模拟退火并不一定能产生绝对的最优解。然而，我们可以按照我们的希望控制其停在当前最小值（和绝对的最优解不同），同时与暴力相比，计算时间可得到很大的缩减。

虚设点有31个可能位置作为初始条件。该模拟退火程序的输入是已知点位置、steiner点的可能位置，还有一个初始路径。（我们手动找了一个最优解，并将其作为程序的输入。我们把已知网络按照三点一组的规则进行分组，每组添加一个虚设点使该组最优，然后使各组最优，得到的配置作为模拟退火的最小路径出现了几次。）

现在配置改变的方式随机从下边选择：

  • 随机在可能位置添加一个steiner点
  • 移除一个已加入的虚设点
  • 将一个已加入的steiner点随机移到一个新的位置

（通过这三种方式，我们可以尝试所有可能的网络配置。如果是为了这样做，第三种方式并没有必要；但是为了给模拟退火步骤更大的自由性，我们包含了这种方法。）

使用存储起来的新配置，为了建立最小生成树和评价网络的花费，一个计算程序会被调用。这个花费然后会被程序的模拟退火部分（routine METROP from Press et al.）使用，以决定新的配置是否被保存或拒绝（根据一个意见一致的冷却计划）。

刚开始的前几次表明一个趋势：增加度为2的虚设点是没用的。虽然虚设点并不花费什么，我们决定在程序中排除这个多余，通过引入较小的虚设点花费，这个花费足够小，并不会影响在必要时加入虚设点，反而最终会处理没有必要的虚设点。

我们同时也添加了一个模拟退火中不经常出现的特点，借此这个程序可以存储测试出的最优解。（模拟退火一般只返回最终使用的配置）

结果

平均迭代9800次退火，1.5min（在25-MHz 386-based PC）

用不同的种子形成不同的随机数生成器，给了5种不同的最优解，长度为94。

超过100个回合时，模拟退火总是收敛到5种中的1种，证明它适用于更大范围的steiner树的问题（当暴力是不可能的时候）。

扩展

如果所有站点是有花费的

如果想用模拟退火，我们可以：

  • 用一些暴力方法
  • 模拟退火数量、位置、虚设点的度（嵌套退火）
  • 修改最小生成树算法，找到花费最小（而不是路径长度最短）的树，并且用前边的模拟退火

似乎前两种方法太费时间而且效率不高，我们选择第三种。

虽然模拟退火不一定能给出最优解，但它在100次试验之后，确实给了最优解（暴力检查过）。

工作记录

  • 2019.4.4

    • 工作

      尝试做第一题，写了一些假设，思路和需要注意的点等等。

    • 感想

      平常还是要做好充足的准备，不要等着题发下来了才去详细了解某些知识和算法。比如这次校赛，昨天出的题，B题就是课上讲过的东西，如果当时抽时间具体研究了，现在不是就已经会了吗？（虽然现在刚学习建模不久，而且当时好像没时间……当然是选择原谅自己）

  • 2019.4.5

    • 工作

      实现第一题matlab程序，发现作的图好丑，不够清晰。寻求解决办法，了解到网络分析工具箱、NetWorkX等。

    • 感想

      擒贼先擒王。最终目的是写一篇好的论文。现在应该把整个思路给定下来，大体上进行实现。不是能因为这个图花费太多时间，应该在后期优化。

  • 2019.4.6

    • 工作

      实现第一题lingo程序，得到答案。重新回顾第一题，整理思路、材料和之后需要做的相关工作。

      看第二题，最小费用最小流，好像得出答案比较简单。了解steiner树，查阅文献。

    • 感想

      道阻且长，我应该对lingo、matlab极其熟悉吗？

  • 2019.4.7

    • 工作

      阅读英文论文点评并总结，着手查看比较好的一篇论文。

    • 感想

      类似的问题，可以直接去看论文。另一方面，更好的是去看论文点评，能够站在顶层往下学习还是很不错的一种体验。

      问题让求steiner树，我们可能是要先证明每个图中steiner树是一定存在的？

  • 2019.4.8

    • 工作

      阅读英文论文模拟退火部分并总结，解决第二题最小费用最大流。

  • 2019.4.9

    • 工作

      上完建模课和队友讨论，知道了凸包..…突然发现题目要求很变态的一点，管道必须是平行于x轴或y轴的。这个是求出结果后不证自明，还是过程中就要求判断呢。

  • 2019.4.10

    • 工作

      知道凸包是什么了，知道了美赛91为什么用L距离，又了解了几种求近似解的方法。着手模拟退火算法

    • 感想

      啊..……不会………..…..…不想说话

  • 2019.4.11

    • 工作

      学习模拟退火，学了一半吧

  • 2019.4.12

    • 工作

      继续学习模拟退火，学完了，知道是什么，大概知道怎么用。

      着手将题目与模拟退火结合，设计代码。

  • 2019.4.13

    • 工作

      敲代码，整体结构敲出来了，差一个扰动函数

    • 感想

      matlab的基础还是不足，很多东西都需要搜着做

  • 2019.4.14

    • 工作

      写扰动函数，尝试调参。难呐..…

  • 2019.4.17

    • 工作

      • 通过暴力，发现这个题和普通的steiner树不一样，按普通的steiner树求解的话，是不符合题目要求的。

        核心呐：它要求边水平或竖直，这样目标函数就应该改成合法性了。

      • 设想了目标函数为合法性的算法如何实现。

  • 2019.4.20

    • 工作

      实现目标函数为合法性的程序，最终可以求得合法的steiner树。但无用的steiner点太多了，需要优化。

  • 2019.4.21

    • 工作

      在合法树的基础上进行长度和steiner点个数的优化。

      遗留了一个问题：自己写的一个函数总是出问题，不符预期。

  • 2019.4.22

    • 工作

      解决昨天遗留的问题，写函数优化无用的steiner点。（原因是在for循环中数组变小了，然后就越界了，python中也有这样的东西嘛我记得

      求得最终答案（不知道强迫症还是什么，总感觉有隐患..…，不过本来的思路就是求近似最优解..）

      整理代码，完善注释。

      那个无用的steiner点即使写了优化函数之后还是有，原因是去掉无用点之后，调用最小生成树算法，这个算法导致有多余的steiner点的产生。

  • 2019.4.23

    • 工作

      讨论结果和论文思路，分工

  • 2019.4.24

    • 工作

      中午下课吃完饭回来就写论文，晚上剪头，吃麻辣烫看了两集BigBang，回来继续写，写到现在11点半了，写了2500+字。主体部分已经写了有四分之三了，一些细节写在了注释里。耶

  • 2019.4.25

    • 工作

      上下午第四节都在写论文，算法设计与分析课上也在写论文..…整体上是写完了，剩下的工作是补充细节。

  • 2019.4.26

    • 工作

      快中午的时候开始补充论文细节，一直到晚上，效率挺低的，可能是分工问题，另外是还有一些讨论的工作。

  • 2019.4.27

    • 工作

      流程图、公式、标题、字体、程序代码等等都搞定了，差不多结束了，第一题代码太烂就先算了吧，那个不重要。

  • 2019.4.28

    • 工作

      搞细节咯，还修改了个流程图。检查了几遍，结束咯~

整个过程想来还是挺艰辛的，大多都是在尝试，不断地试错然后改进，可能建模都这样吧。

作者：@臭咸鱼

转载请注明出处：https://chouxianyu.github.io

欢迎讨论和交流!