数学建模：图论模型 — 最短路模型示例 (Python 求解)

例 1: 设备更新问题

某种工程设备的役龄为 4 年, 每年年初都面临着是否更新的问题: 若卖旧买新, 就要支付一定的购置费用; 若继续使用, 则要支付更多的维护费用, 且使用年限越长维护费用越多. 役龄期内每年的年初购置价格, 当年维护费用及年末剩余净值如下表所示. 为该设备制定一个 4 年役龄期内的更新计划, 使总的支付费用最少.

年份	1	2	3	4
年初购置价格 (万元）	25	26	28	31
当年维护费用 (万元）	10	14	18	26
年末剩余净值 (万元）	20	16	13	11

可以把这个问题化为图论中的最短路问题.

构造赋权有向图 $D=(V,A,W)$, 其中顶点集 V = { v 1 , v 2 , ⋯   , v 5 } V=\left\{v_{1}, v_{2}, \cdots, v_{5}\right\} V={v1​,v2​,⋯,v5​}, 这里 $v_i$ $(i=1,2,3,4)$ 表示第 $i$ 年年初的时刻, $v_5$ 表示第 4 年年末的时刻, $A$ 为边集, 邻接矩阵 $W={\left(w_{ij}\right)}_{5\times 5}$, 这里 $w_{ij}$ 为第 $i$ 年年初至第 $j$ 年年初 (或 $j-1$ 年年末) 期间所支付的费用, 计算公式为
$$w_{ij}=p_i+\sum _{k=1}^{j-i}{a}_k-r_{j-i},$$

其中 $p_i$ 为第 $i$ 年年初的购置价格, $a_k$ 为使用到第 $k$ 年当年的维护费用, $r_i$ 为使用 $i$ 年旧设备的出售价格.

则邻接矩阵
$$W=\left[\begin{array}{ccccc}0& 15& 33& 54& 82\\ \infty & 0& 16& 34& 55\\ \infty & \infty & 0& 18& 36\\ \infty & \infty & \infty & 0& 21\\ \infty & \infty & \infty & \infty & 0\end{array}\right].$$

则制定总的支付费用最小的设备更新计划, 就是求有向图 $D$ 从顶点 $v_1$ 到页点 $v_5$ 的最短路.

import numpy as np
import networkx as nx
import pylab as plt
p=[25,26,28,31]; a=[10,14,18,26]; r=[20,16,13,11]
b=np.zeros((5,5)); #初始化邻接矩阵
for i in range(5):
    for j in range(i+1,5):
        b[i,j]=p[i]+np.sum(a[0:j-i])-r[j-i-1];
print(b)
G=nx.DiGraph(b)
p=nx.dijkstra_path(G, source=0, target=4, weight='weight')  #求最短路径；
print("最短路径为:",np.array(p)+1)  #python下标从0开始
d=nx.dijkstra_path_length(G, 0, 4, weight='weight') #求最短距离
print("所求的费用最小值为：",d)
s=dict(zip(range(5),range(1,6))) #构造用于顶点标注的标号字典
plt.rc('font',size=16)
pos=nx.shell_layout(G)  #设置布局
w=nx.get_edge_attributes(G,'weight')
nx.draw(G,pos,font_weight='bold',labels=s,node_color='r')
nx.draw_networkx_edge_labels(G,pos,edge_labels=w)
path_edges=list(zip(p,p[1:]))
nx.draw_networkx_edges(G,pos,edgelist=path_edges,edge_color='r',width=3) #画出最短路径
#plt.savefig("figure10_9.png",pdi=500); 
plt.show()

[[ 0. 15. 33. 54. 82.]
 [ 0.  0. 16. 34. 55.]
 [ 0.  0.  0. 18. 36.]
 [ 0.  0.  0.  0. 21.]
 [ 0.  0.  0.  0.  0.]]
最短路径为: [1 2 3 5]
所求的费用最小值为： 67.0

求得 $v_1$ 到 $v_5$ 的最短路径为 $v_1\to v_2\to v_3\to v_5$, 最短路径的长度为67. 即设备更新计划为第1年年初买进新设备, 使用到第 1 年年底, 第 2 年年初购进新设备, 使用到第 2 年年底, 第 3 年年初再购进新设备, 使用到第 4 年年底.

例 2: 重心问题

重心问题指有些公共服务设施 (例如邮局, 学校等) 的选址, 要求设施到所有服务对象点的距离总和最小. 一般要考虑人口密度问题, 或者全体被服务对象来往的总路程最短. 例如下面的问题:

某矿区有六个产矿点, 如图所示, 已知各产矿点每天的产矿量 (标在图中的各顶点旁) 为 $q_i(i=1,2,\cdots ,6)$ 吨, 现要从这六个产矿点选一个来建造选矿厂, 问应选在哪个产矿点, 才能使各产矿点所产的矿石运到选矿厂所在地的总运力 (吨·公里) 最小.

令 $d_{ij}(i,j=1,2,\cdots ,6)$ 表示顶点 $v_i$ 与 $v_j$ 之间的距离. 若选矿厂设在 $v_i$ 并且各产矿点到选矿厂的总运力为 $m_i$, 则确定选矿厂的位置就转化为求 $m_k$, 使得 $m_k={\min }_{1\le i\le 6}{m}_i$.

由于各产矿点到选矿厂的总运力依赖于任意两顶点之间的距离, 即任意两顶点之间最短路的长度, 因此可首先利用 Dijkstra (或 Floyd) 算法求出所有顶点对之间的最短距离, 然后计算出顶点 $v_i$ 设立选矿厂时各产矿点到 $v_i$ 的总运力
$$m_i=\sum _{j=1}^6{q}_j{d}_{ij},i=1,2,\cdots ,6.$$

最后利用 ${\min }_{1\le i\le 6}{m}_i$ 选择设置选矿厂的位置.

计算的 Python 程序如下:

import numpy as np
import networkx as nx
List=[(1,2,20),(1,5,15),(2,3,20),(2,4,40),
      (2,5,25),(3,4,30),(3,5,10),(5,6,15)]
G=nx.Graph()
G.add_nodes_from(range(1,7))
G.add_weighted_edges_from(List)
c=dict(nx.shortest_path_length(G,weight='weight'))
d=np.zeros((6,6))
for i in range(1,7):
    for j in range(1,7): d[i-1,j-1]=c[i][j]
print(d)
q=np.array([80,90,30,20,60,10])
m=d@q  #计算运力，这里使用矩阵乘法
mm=m.min()  #求运力的最小值
ind=np.where(m==mm)[0]+1  #python下标从0开始，np.where返回值为元组
print("运力m=",m,'\n最小运力mm=',mm,"\n选矿厂的设置位置为：",ind)

[[ 0. 20. 25. 55. 15. 30.]
 [20.  0. 20. 40. 25. 40.]
 [25. 20.  0. 30. 10. 25.]
 [55. 40. 30.  0. 40. 55.]
 [15. 25. 10. 40.  0. 15.]
 [30. 40. 25. 55. 15.  0.]]
运力m= [ 4850.  4900.  5250. 11850.  4700.  8750.] 
最小运力mm= 4700.0 
选矿厂的设置位置为： [5]

${\min }_{1\le i\le 6}{m}_i=m_5$, 所以 $v_5$ 为设置选矿厂的位置