司守奎《数学建模算法与应用》第二版 第四章图论习题答案
4.1
import networkx as nx
import matplotlib.pyplot as plt
if __name__ == "__main__":
graph = nx.Graph() # 创建无向图
graph.add_edge("L", "M", weight=56)
graph.add_edge("L", "N", weight=35)
graph.add_edge("L", "Pa", weight=21)
graph.add_edge("L", "Pe", weight=51)
graph.add_edge("L", "T", weight=60)
graph.add_edge("M", "N", weight=21)
graph.add_edge("M", "Pa", weight=57)
graph.add_edge("M", "Pe", weight=78)
graph.add_edge("M", "T", weight=70)
graph.add_edge("N", "Pa", weight=36)
graph.add_edge("N", "Pe", weight=68)
graph.add_edge("N", "T", weight=68)
graph.add_edge("Pa", "Pe", weight=51)
graph.add_edge("Pa", "T", weight=61)
graph.add_edge("Pe", "T", weight=13)
res = nx.minimum_spanning_tree(graph) # 使用Kruskal算法计算最小生成树
print("最小生成树为",res.edges)
print("最小生成树的权重和为", sum(d["weight"] for u,v,d in res.edges(data=True)))
# 画图
plt.figure(figsize=(8, 6)) # 新建画布
pos = nx.spring_layout(graph)
nx.draw(graph, pos, with_labels=True, node_size=500, node_color="skyblue", font_size=12, font_weight="bold") # 绘制无向图的节点
nx.draw_networkx_edge_labels(graph, pos, edge_labels={(u, v): d["weight"] for u, v, d in graph.edges(data=True)}) # 绘制无向图的边界
nx.draw_networkx_edges(res, pos, edge_color="red", width=2)# 绘制最小生成树的边界
plt.title("Minimum Spanning Tree")
plt.show()
# output:最小生成树为 [('L', 'Pa'), ('L', 'N'), ('L', 'Pe'), ('M', 'N'), ('Pe', 'T')]
# 最小生成树的权重和为 141
绘制出的图形如下
4.2
假设
v
i
,
i
=
1
,
2
,
3
,
4
,
5
v_i,i=1,2,3,4,5
vi,i=1,2,3,4,5表示第
i
i
i年初。
w
i
j
w_{ij}
wij表示第
i
i
i年到第
j
j
j年的费用（机器第
i
i
i年购入，第
j
j
j年卖出）。
w
12
=
0.3
+
2.5
−
2.0
=
0.8
w
13
=
0.3
+
0.8
+
2.5
−
1.6
=
2.0
w
14
=
0.3
+
0.8
+
1.5
+
2.5
−
1.3
=
3.8
w
15
=
0.3
+
0.8
+
1.5
+
2.0
+
2.5
−
1.1
=
6
w
23
=
0.3
+
2.6
−
2.0
=
0.9
w
24
=
0.3
+
0.8
+
2.6
−
1.6
=
2.1
w
25
=
0.3
+
0.8
+
1.5
+
2.6
−
1.3
=
3.9
w
34
=
0.3
+
2.8
−
2.0
=
1.1
w
35
=
0.3
+
0.8
+
2.8
−
1.6
=
2.3
w
45
=
0.3
+
3.1
−
2.0
=
1.4
w_{12}=0.3+2.5-2.0=0.8\\ w_{13}=0.3+0.8+2.5-1.6=2.0\\ w_{14}=0.3+0.8+1.5+2.5-1.3=3.8\\ w_{15}=0.3+0.8+1.5+2.0+2.5-1.1=6\\ w_{23}=0.3+2.6-2.0=0.9\\ w_{24}=0.3+0.8+2.6-1.6=2.1\\ w_{25}=0.3+0.8+1.5+2.6-1.3=3.9\\ w_{34}=0.3+2.8-2.0=1.1\\ w_{35}=0.3+0.8+2.8-1.6=2.3\\ w_{45}=0.3+3.1-2.0=1.4
w12=0.3+2.5−2.0=0.8w13=0.3+0.8+2.5−1.6=2.0w14=0.3+0.8+1.5+2.5−1.3=3.8w15=0.3+0.8+1.5+2.0+2.5−1.1=6w23=0.3+2.6−2.0=0.9w24=0.3+0.8+2.6−1.6=2.1w25=0.3+0.8+1.5+2.6−1.3=3.9w34=0.3+2.8−2.0=1.1w35=0.3+0.8+2.8−1.6=2.3w45=0.3+3.1−2.0=1.4
实际上，问题也就是从节点 v 1 v_1 v1到节点 v 5 v_5 v5的最小距离。
import networkx as nx
import matplotlib.pyplot as plt
if __name__ == "__main__":
graph = nx.DiGraph() # 创建有向图
graph.add_edge("v1", "v2", weight=0.8)
graph.add_edge("v1", "v3", weight=2)
graph.add_edge("v1", "v4", weight=3.8)
graph.add_edge("v1", "v5", weight=6)
graph.add_edge("v2", "v3", weight=0.9)
graph.add_edge("v2", "v4", weight=2.1)
graph.add_edge("v2", "v5", weight=3.9)
graph.add_edge("v3", "v4", weight=1.1)
graph.add_edge("v3", "v5", weight=2.3)
graph.add_edge("v4", "v5", weight=1.4)
res = nx.shortest_path(
graph, "v1", "v5", weight="weight", method="dijkstra"
) # 使用dijkstra算法求解最短路径
# 计算指定路径的权重和
total_weight = 0
for i in range(len(res) - 1):
u = res[i]
v = res[i + 1]
total_weight += graph[u][v]["weight"]
print("最短路径为", res)
print("最短路径长度为", total_weight)
# 画图
plt.figure(figsize=(8, 6))
pos = nx.spring_layout(graph)
nx.draw(graph, pos, with_labels=True, node_size=500, node_color="skyblue", font_size=12, font_weight="bold") # 绘制有向图的节点
nx.draw_networkx_edge_labels(graph, pos, edge_labels={(u, v): d["weight"] for u, v, d in graph.edges(data=True)}) # 绘制整个有向图的边界
nx.draw_networkx_edges(graph, pos, edgelist=[(res[i], res[i + 1]) for i in range(len(res) - 1)], edge_color="red", width=2) # 绘制最短路径边界
plt.title("Shortest Path")
plt.show()
# output:最短路径为 ['v1', 'v2', 'v3', 'v5']
# 最短路径长度为 4.0
绘制的图形如下
4.3
import networkx as nx
import matplotlib.pyplot as plt
if __name__=="__main__":
graph = nx.DiGraph() # 创建有向图
graph.add_edge('s','A',capacity=20)
graph.add_edge('s','B',capacity=20)
graph.add_edge('s','C',capacity=100)
graph.add_edge('A', '1', capacity=30)
graph.add_edge('A','2',capacity=10)
graph.add_edge('A','4',capacity=40)
graph.add_edge('B','3',capacity=10)
graph.add_edge('B','4',capacity=50)
graph.add_edge('C','1',capacity=20)
graph.add_edge('C','2',capacity=10)
graph.add_edge('C','3',capacity=40)
graph.add_edge('C','4',capacity=5)
graph.add_edge('1','t',capacity=20)
graph.add_edge('2','t',capacity=20)
graph.add_edge('3','t',capacity=60)
graph.add_edge('4','t',capacity=20)
res = nx.maximum_flow(graph, 's','t')
print("网络最大流为", res[0])
# 画图
plt.figure(figsize=(8, 6))
pos = nx.spring_layout(graph)
nx.draw(graph, pos, with_labels=True, node_size=500, node_color="skyblue", font_size=12, font_weight="bold") # 绘制有向图的节点
nx.draw_networkx_edge_labels(graph, pos, edge_labels={(u, v): d["capacity"] for u, v, d in graph.edges(data=True)}) # 绘制整个有向图的标签
nx.draw_networkx_edges(graph, pos, edgelist=[(src,dst) for i in range(1,len(res)) for src,j in res[i].items() for dst in j.keys()], edge_color="red", width=2) # 绘制网络最大流边界
nx.draw_networkx_edge_labels(graph, pos,edge_labels={(src,dst):k for i in range(1,len(res)) for src,j in res[i].items() for dst,k in j.items()},label_pos=0.8,font_color="red") # 绘制网络最大流标签
plt.title("Maximum Flow")
plt.show()
# output:网络最大流为 110
绘制的网络最大流图如下
由图可以看出，只有市场3不满足，只能满足50单位，还差10单位。
4.4
假设甲、乙、丙、丁、戊分别为 a , b , c , d , e a,b,c,d,e a,b,c,d,e，俄文、英文、日文、德文、法文分别为 1 , 2 , 3 , 4 , 5 1,2,3,4,5 1,2,3,4,5，则
import networkx as nx
import matplotlib.pyplot as plt
if __name__ == "__main__":
graph = nx.DiGraph()
graph.add_edge("s", "a", capacity=1)
graph.add_edge("s", "b", capacity=1)
graph.add_edge("s", "c", capacity=1)
graph.add_edge("s", "d", capacity=1)
graph.add_edge("s", "e", capacity=1)
graph.add_edge("a", "2", capacity=1)
graph.add_edge("a", "3", capacity=1)
graph.add_edge("b", "1", capacity=1)
graph.add_edge("b", "2", capacity=1)
graph.add_edge("b", "4", capacity=1)
graph.add_edge("c", "2", capacity=1)
graph.add_edge("c", "3", capacity=1)
graph.add_edge("d", "3", capacity=1)
graph.add_edge("d", "2", capacity=1)
graph.add_edge("e", "4", capacity=1)
graph.add_edge("e", "5", capacity=1)
graph.add_edge("1", "t", capacity=1)
graph.add_edge("2", "t", capacity=1)
graph.add_edge("3", "t", capacity=1)
graph.add_edge("4", "t", capacity=1)
graph.add_edge("5", "t", capacity=1)
res = nx.maximum_flow(graph, "s", "t")
print("网络最大流为", res[0])
# 画图
plt.figure(figsize=(8, 6))
pos = nx.spring_layout(graph)
nx.draw(
graph,
pos,
with_labels=True,
node_size=500,
node_color="skyblue",
font_size=12,
font_weight="bold",
) # 绘制有向图的节点
nx.draw_networkx_edge_labels(
graph,
pos,
edge_labels={(u, v): d["capacity"] for u, v, d in graph.edges(data=True)},
) # 绘制整个有向图的标签
nx.draw_networkx_edges(
graph,
pos,
edgelist=[
(src, dst)
for i in range(1, len(res))
for src, j in res[i].items()
for dst in j.keys()
],
edge_color="red",
width=2,
) # 绘制网络最大流边界
nx.draw_networkx_edge_labels(
graph,
pos,
edge_labels={
(src, dst): k
for i in range(1, len(res))
for src, j in res[i].items()
for dst, k in j.items()
},
label_pos=0.8,
font_color="red",
) # 绘制网络最大流标签
plt.title("Maximum Flow")
plt.show()
# output:网络最大流为 4
绘制的图形如下
由上图可知，从 s s s点出发，只有 a a a也即是甲的流量为0，所以只有甲没有收到聘书。乙获得了俄语的聘书，丙获得了英语的聘书，丁获得了日语的聘书，戊获得了德语的聘书。
4.5
import networkx as nx
import matplotlib.pyplot as plt
if __name__=="__main__":
graph = nx.DiGraph()
graph.add_edge('s','a',capacity=8,weight=0)
graph.add_edge('s','b',capacity=7,weight=0)
graph.add_edge('a','1',capacity=8,weight=20)
graph.add_edge('a','2',capacity=8,weight=24)
graph.add_edge('a','3',capacity=8,weight=5)
graph.add_edge('b','1',capacity=7,weight=30)
graph.add_edge('b','2',capacity=7,weight=22)
graph.add_edge('b','3',capacity=7,weight=20)
graph.add_edge('1','t',capacity=4,weight=0)
graph.add_edge('2','t',capacity=5,weight=0)
graph.add_edge('3','t',capacity=6,weight=0)
res = nx.max_flow_min_cost(graph, 's','t')
mincost = nx.cost_of_flow(graph, res)
print("网络流的最小花费为", mincost)
# 画图
plt.figure(figsize=(8, 6))
pos = nx.spring_layout(graph)
nx.draw(
graph,
pos,
with_labels=True,
node_size=500,
node_color="skyblue",
font_size=12,
font_weight="bold",
) # 绘制有向图的节点
nx.draw_networkx_edge_labels(
graph,
pos,
edge_labels={(u, v): (d["capacity"],d["weight"]) for u, v, d in graph.edges(data=True)},
) # 绘制整个有向图的标签
nx.draw_networkx_edges(
graph,
pos,
edgelist=[
(src, dst)
for src,i in res.items()
for dst in i.keys()
],
edge_color="red",
width=2,
) # 绘制网络最大流边界
nx.draw_networkx_edge_labels(
graph,
pos,
edge_labels={
(src, dst): k
for src,i in res.items()
for dst,k in i.items()
},
label_pos=0.8,
font_color="red",
) # 绘制网络最大流标签
plt.title("Maximum Flow")
plt.show()
# output:网络流的最小花费为 240
绘制的网络流图如下
4.6
import networkx as nx
import matplotlib.pyplot as plt
if __name__=="__main__":
graph = nx.DiGraph()
graph.add_edge('vs', 'v1',capacity=10,weight=4)
graph.add_edge('vs','v2',capacity=8,weight=1)
graph.add_edge('v1','vt',capacity=7,weight=1)
graph.add_edge('v2','v1',capacity=5,weight=2)
graph.add_edge('v2','v3',capacity=10,weight=3)
graph.add_edge('v1','v3',capacity=2,weight=6)
graph.add_edge('v3','vt',capacity=4,weight=2)
res=nx.max_flow_min_cost(graph,'vs','vt')
mincost = nx.cost_of_flow(graph, res)
print("网络流的最小花费为", mincost)
# 画图
plt.figure(figsize=(8, 6))
pos = nx.spring_layout(graph)
nx.draw(
graph,
pos,
with_labels=True,
node_size=500,
node_color="skyblue",
font_size=12,
font_weight="bold",
) # 绘制有向图的节点
nx.draw_networkx_edge_labels(
graph,
pos,
edge_labels={(u, v): (d["capacity"],d["weight"]) for u, v, d in graph.edges(data=True)},
) # 绘制整个有向图的标签
nx.draw_networkx_edges(
graph,
pos,
edgelist=[
(src, dst)
for src,i in res.items()
for dst in i.keys()
],
edge_color="red",
width=2,
) # 绘制网络最大流边界
nx.draw_networkx_edge_labels(
graph,
pos,
edge_labels={
(src, dst): k
for src,i in res.items()
for dst,k in i.items()
},
label_pos=0.8,
font_color="red",
) # 绘制网络最大流标签
plt.title("Maximum Flow")
plt.show()
# output:网络流的最小花费为 55
绘制的图形如下
4.7
（1）
import networkx as nx
import matplotlib.pyplot as plt
if __name__ == "__main__":
# 创建一个有向图
G = nx.DiGraph()
# 添加活动节点
activities = {
"A": {"duration": 6},
"B": {"duration": 5},
"C": {"duration": 3},
"D": {"duration": 2},
"E": {"duration": 3},
"F": {"duration": 2},
"G": {"duration": 4},
"H": {"duration": 2},
}
for activity, attr in activities.items():
G.add_node(activity, **attr)
# 添加活动之间的依赖关系
dependencies = [
("A", "C"),
("C", "D"),
("A", "E"),
("D", "E"),
("B", "F"),
("E", "G"),
("G", "H"),
("F", "H"),
]
G.add_edges_from(dependencies)
# 绘制网络图
pos = nx.spring_layout(G)
nx.draw(
G,
pos,
with_labels=True,
node_size=700,
node_color="skyblue",
font_size=12,
font_weight="bold",
arrowsize=20,
)
plt.title("Project Network Diagram")
plt.show()
（2）
求最早开工时间 E S ES ES的时候，不同于教材里面使用线性规划的方法，这里采用了类似与贪心算法的求解。关键路径这里有一个误区是，并不是说完成这个工程只需要进行关键路径上的事件就完事了，而是所有的事件都要进行，但是关键路径是耗时最长的事件线，所以当其他事件线在并行执行的时候，我们只需要考虑关键路径的耗时来进行规划，这也解释了为什么线性规划中需要使用 min \min min进行求解，而在这里使用了 max \max max进行求解。我们要求的最小，是在完成所有事件的前提下，最长的那条关键路径的最小。
假设
x
i
x_i
xi是事件
i
i
i开始的最早时间，
z
i
z_i
zi是开始的最晚时间，
t
i
j
t_{ij}
tij是作业
(
i
,
j
)
(i,j)
(i,j)的计划时间，
e
s
i
j
,
l
s
i
j
,
e
f
i
j
,
l
f
i
j
es_{ij},ls_{ij},ef_{ij},lf_{ij}
esij,lsij,efij,lfij分别表示作业
(
i
j
)
(ij)
(ij)的最早开始时间，最迟开工时间，最早完工时间，最晚完工时间。因此有
z
n
=
x
n
z
i
=
min
j
{
z
j
−
t
i
j
}
,
i
=
n
−
1
,
.
.
,
1
,
(
i
,
j
)
∈
A
e
s
i
j
=
x
i
,
(
i
,
j
)
∈
A
l
f
i
j
=
z
j
,
(
i
,
j
)
∈
A
l
s
i
j
=
l
f
i
j
−
t
i
j
,
(
i
,
j
)
∈
A
e
f
i
j
=
x
i
+
t
i
j
,
(
i
,
j
)
∈
A
z_n=x_n\\ z_i=\min\limits_{j}\{z_j-t_{ij}\},i=n-1,..,1,(i,j)\in A\\ es_{ij}=x_i,(i,j)\in A\\ lf_{ij}=z_j,(i,j)\in A\\ ls_{ij}=lf_{ij}-t_{ij},(i,j)\in A\\ ef_{ij}=x_i+t_{ij},(i,j)\in A
zn=xnzi=jmin{zj−tij},i=n−1,..,1,(i,j)∈Aesij=xi,(i,j)∈Alfij=zj,(i,j)∈Alsij=lfij−tij,(i,j)∈Aefij=xi+tij,(i,j)∈A
代码如下
import networkx as nx
import matplotlib.pyplot as plt
if __name__ == "__main__":
# 创建一个有向图
G = nx.DiGraph()
# 添加活动节点
activities = {
"A": {"duration": 6},
"B": {"duration": 5},
"C": {"duration": 3},
"D": {"duration": 2},
"E": {"duration": 3},
"F": {"duration": 2},
"G": {"duration": 4},
"H": {"duration": 2},
}
for activity, attr in activities.items():
G.add_node(activity, **attr)
# 添加活动之间的依赖关系
dependencies = [
("A", "C"),
("C", "D"),
("A", "E"),
("D", "E"),
("B", "F"),
("E", "G"),
("G", "H"),
("F", "H"),
]
G.add_edges_from(dependencies)
# 计算最早开始时间和最早完成时间
es = list(nx.topological_sort(G)) # 返回的是一个迭代器，转换成列表形式以供后续使用
ef = {node: 0 for node in es}
for node in es:
max_ef = 0
if (
len(list(G.predecessors(node))) != 0
): # 如果当前node的前驱节点存在，则更新，否则不更新
for predecessor in G.predecessors(node):
max_ef = max(max_ef, ef[predecessor] + G.nodes[predecessor]["duration"])
ef[node] = max_ef
# 计算最迟开始时间和最迟完成时间
lf = {node: ef[node] for node in ef}
for node in reversed(list(nx.topological_sort(G))):
min_lf = float("inf")
if (
len(list(G.successors(node))) != 0
): # 如果当前node的后驱节点存在，则更新，否则不更新
for successor in G.successors(node):
min_lf = min(min_lf, lf[successor] - G.nodes[node]["duration"])
lf[node] = min_lf
# 计算活动的浮动时间
float_time = {node: lf[node] - ef[node] for node in ef}
# 计算关键路径
critical_path = [
(u, v) for u, v in G.edges() if lf[v] - ef[u] == G.nodes[u]["duration"]
]
# 打印结果
for node in sorted(activities.keys()):
print(
f"事件{node}的最早开始时间{ef[node]}最早结束时间{ef[node] + activities[node]['duration']}最晚开始时间{lf[node]}最晚结束时间{lf[node] + activities[node]['duration']}"
)
print("关键路径为", critical_path)
# 绘制网络图
pos = nx.spring_layout(G)
nx.draw(
G,
pos,
with_labels=True,
node_size=700,
node_color="skyblue",
font_size=12,
font_weight="bold",
arrowsize=20,
)
plt.title("Project Network Diagram")
plt.show()
# output:事件A的最早开始时间0最早结束时间6最晚开始时间0最晚结束时间6
# 事件B的最早开始时间0最早结束时间5最晚开始时间11最晚结束时间16
# 事件C的最早开始时间6最早结束时间9最晚开始时间6最晚结束时间9
# 事件D的最早开始时间9最早结束时间11最晚开始时间9最晚结束时间11
# 事件E的最早开始时间11最早结束时间14最晚开始时间11最晚结束时间14
# 事件F的最早开始时间5最早结束时间7最晚开始时间16最晚结束时间18
# 事件G的最早开始时间14最早结束时间18最晚开始时间14最晚结束时间18
# 事件H的最早开始时间18最早结束时间20最晚开始时间18最晚结束时间20
# 关键路径为 [('A', 'C'), ('C', 'D'), ('D', 'E'), ('E', 'G'), ('G', 'H')]
绘制的关键路径如下
（3）
假设
x
i
x_i
xi是事件
i
i
i的开始时间，
t
i
j
t_{ij}
tij是作业
(
i
,
j
)
(i,j)
(i,j)的计划时间，
m
i
j
m_{ij}
mij是完成作业
(
i
,
j
)
(i,j)
(i,j)的最短时间，
y
i
j
y_{ij}
yij是作业
(
i
,
j
)
(i,j)
(i,j)可能减少的时间，
c
i
j
c_{ij}
cij是作业
(
i
,
j
)
(i,j)
(i,j)缩短一天增加的费用，限制完成的时间为
d
d
d，因此有
x
j
−
x
i
≥
t
i
j
−
y
i
j
且
0
≤
y
i
j
≤
t
i
j
−
m
i
j
min
∑
(
i
,
j
)
∈
A
c
i
j
y
i
j
s
.
t
.
{
x
j
−
x
i
+
y
i
j
≥
t
i
j
,
(
i
,
j
)
∈
A
,
i
,
j
∈
V
x
n
−
x
1
≤
d
0
≤
y
i
j
≤
t
i
j
−
m
i
j
,
(
i
,
j
)
∈
A
,
i
,
j
∈
V
x_j-x_i\geq t_{ij}-y_{ij}且0\leq y_{ij}\leq t_{ij}-m_{ij}\\ \min\sum\limits_{(i,j)\in A}c_{ij}y_{ij}\\ \begin{equation} s.t. \begin{cases} x_j-x_i+y_{ij}\geq t_{ij},(i,j)\in A,i,j\in V\\ x_n-x_1\leq d\\ 0\leq y_{ij}\leq t_{ij}-m_{ij},(i,j)\in A,i,j\in V \end{cases} \end{equation}
xj−xi≥tij−yij且0≤yij≤tij−mijmin(i,j)∈A∑cijyijs.t.⎩
⎨
⎧xj−xi+yij≥tij,(i,j)∈A,i,j∈Vxn−x1≤d0≤yij≤tij−mij,(i,j)∈A,i,j∈V
代码如下
import networkx as nx
import numpy as np
from scipy.optimize import linprog
if __name__ == "__main__":
f=np.zeros(16) # 前8位表示最早开始时间
f[8:]=[800,600,300,600,400,300,200,200] # 后8位表示减少的周数
A=np.array([[1,-1,0,0,0,0,0,0,-1,0,0,0,0,0,0,0],
[1,0,-1,0,0,0,0,0,0,-1,0,0,0,0,0,0],
[0,1,0,-1,0,0,0,0,0,0,-1,0,0,0,0,0],
[0,0,0,1,-1,0,0,0,0,0,0,-1,0,0,0,0],
[0,0,0,0,1,-1,0,0,0,0,0,0,-1,0,0,0],
[0,0,1,0,0,0,-1,0,0,0,0,0,0,-1,0,0],
[0,0,0,0,0,1,-1,0,0,0,0,0,0,0,-1,0],
[0,0,0,0,0,0,1,-1,0,0,0,0,0,0,0,-1],
[0,1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0],
[-1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0]])
b = np.array([-6,-5,-3,-2,-3,-2,-4,-2,0,17]).T
lb = np.zeros(16)
ub = np.array([None]*16)
ub[8:] = [6-4,5-3,3-1,2-1,3-1,2-1,4-2,2-2]
bound = np.vstack((lb,ub)).T
res=linprog(f,A,b,bounds=bound)
print("最优解为",res.x)
print("最优值为",res.fun)
# output最优解为 [ 0. 6. 5. 8. 10. 13. 15. 17. 0. 0. 1. 0. 0. 0. 2. 0.]
# 最优值为 700.0
（4）
假设
t
i
j
,
a
i
j
,
e
i
j
,
b
i
j
t_{ij},a_{ij},e_{ij},b_{ij}
tij,aij,eij,bij分别是完成作业
(
i
,
j
)
(i,j)
(i,j)的实际时间，最乐观时间，最可能时间，最悲观时间，通常用下面的方法计算相应的数学期望和方差：
E
(
t
i
j
)
=
a
i
j
+
4
e
i
j
+
b
i
j
6
v
a
r
(
t
i
j
)
=
(
b
i
j
−
a
i
j
)
2
36
实际工期
T
=
∑
(
i
,
j
)
∈
关键路线
t
i
j
E(t_{ij})=\frac{a_{ij}+4e_{ij}+b_{ij}}{6}\\ var(t_{ij})=\frac{(b_{ij}-a_{ij})^2}{36}\\ 实际工期T=\sum\limits_{(i,j)\in 关键路线}t_{ij}
E(tij)=6aij+4eij+bijvar(tij)=36(bij−aij)2实际工期T=(i,j)∈关键路线∑tij
由中心极限定理，可以假设
T
T
T服从正态分布，并且期望值和方差满足
T
ˉ
=
E
(
T
)
=
∑
(
i
,
j
)
∈
关键路线
E
(
t
i
j
)
S
2
=
v
a
r
(
T
)
=
∑
(
i
,
j
)
∈
关键路线
v
a
r
(
t
i
j
)
\bar{T}=E(T)=\sum\limits_{(i,j)\in 关键路线}E(t_{ij})\\ S^2=var(T)=\sum\limits_{(i,j)\in 关键路线}var(t_{ij})
Tˉ=E(T)=(i,j)∈关键路线∑E(tij)S2=var(T)=(i,j)∈关键路线∑var(tij)
设规定的工期为
d
d
d，则在规定的工期内完成整个项目的概率为
P
{
T
≤
d
}
=
Φ
(
d
−
T
ˉ
S
)
P\{T\leq d\}=\Phi(\frac{d-\bar{T}}{S})
P{T≤d}=Φ(Sd−Tˉ)
使用线性规划来求关键路径
max
∑
(
i
,
j
)
∈
A
E
(
t
i
j
)
x
i
j
s
.
t
.
{
∑
j
:
(
i
,
j
)
∈
A
x
i
j
−
∑
j
(
j
,
i
)
∈
A
x
j
i
=
{
1
,
i
=
1
−
1
,
i
=
n
0
,
i
≠
1
,
n
x
i
j
=
0
或
1
,
(
i
,
j
)
∈
A
\max\sum\limits_{(i,j)\in A}E(t_{ij})x_{ij}\\ \begin{equation} s.t. \begin{cases} \sum\limits_{j:(i,j)\in A}x_{ij}-\sum\limits_{j(j,i)\in A}x_{ji}=\begin{cases}1,&i=1\\-1,&i=n\\0,&i\neq 1,n\end{cases}\\ x_{ij}=0或1,&(i,j)\in A \end{cases} \end{equation}
max(i,j)∈A∑E(tij)xijs.t.⎩
⎨
⎧j:(i,j)∈A∑xij−j(j,i)∈A∑xji=⎩
⎨
⎧1,−1,0,i=1i=ni=1,nxij=0或1,(i,j)∈A
代码如下
import numpy as np
from scipy.stats import norm
from scipy.optimize import linprog
if __name__ == "__main__":
a_ij = np.array([2, 4, 2, 1, 1, 3, 2, 1])
e_ij = np.array([6, 5, 3, 2, 3, 4, 4, 2])
b_ij = np.array([10, 6, 4, 3, 5, 5, 6, 4])
E = (a_ij + 4 * e_ij + b_ij) / 6
var = (b_ij - a_ij) ** 2 / 36
d = 21
p = 0.95
f = -E
Aeq = np.array(
[
[1, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, -1],
[-1, 0, 1, 0, 0, 0, 0, 0],
[0, -1, 0, 0, 0, 1, 0, 0],
[0, 0, -1, 1, 0, 0, 0, 0],
[0, 0, 0, -1, 1, 0, 0, 0],
[0, 0, 0, 0, -1, 0, 1, 0],
[0, 0, 0, 0, 0, -1, -1, 1],
]
)
beq = np.array([[1, -1, 0, 0, 0, 0, 0, 0]]).T
lb = np.zeros(8)
ub = np.ones(8)
bound = np.vstack((lb, ub)).T
res = linprog(f, A_eq=Aeq, b_eq=beq, bounds=bound, integrality=3)
E = sum(E[res.x == 1]) # 关键路线的数学期望
S = sum(var[res.x == 1]) ** (1 / 2) # 关键路线的标准差
cdf = norm.cdf(d, E, S) # 累积分布函数求样本值的概率
ppf = norm.ppf(p, E, S) # 逆累积分布函数求给定概率的样本值
print("在21周上市的概率为{}%".format(round(cdf * 100, 2)))
print("以95%的概率完成新品上市所需要的周数为{}周".format(round(ppf, 2)))
# output:在21周上市的概率为68.04%
# 以95%的概率完成新品上市所需要的周数为23.08周
4.8
假设
v
i
,
i
=
1
,
2
,
3
,
4
,
5
,
6
v_i,i=1,2,3,4,5,6
vi,i=1,2,3,4,5,6表示第
i
i
i年初，
w
i
j
w_{ij}
wij表示第
i
i
i年到第
j
j
j年支付的费用（新机器第
i
i
i年买入，第
j
j
j年卖出）。
w
12
=
11
+
5
=
16
w
13
=
11
+
5
+
6
=
22
w
14
=
11
+
5
+
6
+
8
=
30
w
15
=
11
+
5
+
6
+
8
+
11
=
41
w
16
=
11
+
5
+
6
+
8
+
11
+
20
=
61
w
23
=
11
+
5
=
16
w
24
=
11
+
5
+
6
=
22
w
25
=
11
+
5
+
6
+
8
=
30
w
26
=
11
+
5
+
6
+
8
+
11
=
41
w
34
=
12
+
5
=
17
w
35
=
12
+
5
+
6
=
23
w
36
=
12
+
5
+
6
+
8
=
31
w
45
=
12
+
5
=
17
w
46
=
12
+
5
+
6
=
23
w
56
=
13
+
5
=
18
w_{12}=11+5=16\\ w_{13}=11+5+6=22\\ w_{14}=11+5+6+8=30\\ w_{15}=11+5+6+8+11=41\\ w_{16}=11+5+6+8+11+20=61\\ w_{23}=11+5=16\\ w_{24}=11+5+6=22\\ w_{25}=11+5+6+8=30\\ w_{26}=11+5+6+8+11=41\\ w_{34}=12+5=17\\ w_{35}=12+5+6=23\\ w_{36}=12+5+6+8=31\\ w_{45}=12+5=17\\ w_{46}=12+5+6=23\\ w_{56}=13+5=18\\
w12=11+5=16w13=11+5+6=22w14=11+5+6+8=30w15=11+5+6+8+11=41w16=11+5+6+8+11+20=61w23=11+5=16w24=11+5+6=22w25=11+5+6+8=30w26=11+5+6+8+11=41w34=12+5=17w35=12+5+6=23w36=12+5+6+8=31w45=12+5=17w46=12+5+6=23w56=13+5=18
5年之内所支付的费用最少，也就是求该有向图的最短路径
import networkx as nx
if __name__=="__main__":
graph = nx.DiGraph()
graph.add_edge('1','2',weight=16)
graph.add_edge('1','3',weight=22)
graph.add_edge('1','4',weight=30)
graph.add_edge('1','5',weight=41)
graph.add_edge('1','6',weight=61)
graph.add_edge('2','3',weight=16)
graph.add_edge('2','4',weight=22)
graph.add_edge('2','5',weight=30)
graph.add_edge('2','6',weight=41)
graph.add_edge('3','4',weight=17)
graph.add_edge('3','5',weight=23)
graph.add_edge('3','6',weight=31)
graph.add_edge('4','5',weight=17)
graph.add_edge('4','6',weight=23)
graph.add_edge('5','6',weight=18)
res = nx.shortest_path(graph,'1','6',weight="weight")
total_weight = 0
for i in range(len(res) - 1):
u = res[i]
v = res[i + 1]
total_weight += graph[u][v]["weight"]
print("最短路径为", res)
print("最短路径长度为", total_weight)
# output:最短路径为 ['1', '3', '6']
# 最短路径长度为 53
4.9
假设
x
i
,
i
=
1
,
2
,
3
,
4
,
5
,
6
,
7
,
8
x_i,i=1,2,3,4,5,6,7,8
xi,i=1,2,3,4,5,6,7,8代表相应的工作代号的开始时间，
t
i
j
t_{ij}
tij为正常工作从工作代号
i
i
i到工作代号
j
j
j所需要的时间，
w
i
j
w_{ij}
wij为特急时间，
d
i
j
d_{ij}
dij为减少的时间，
q
i
j
q_{ij}
qij为正常工作时间所花费的费用，
p
i
j
p_{ij}
pij为特急时间所花费的费用。
min
{
15
(
x
8
−
x
1
)
+
∑
(
i
,
j
)
∈
A
q
i
j
(
t
i
j
−
d
i
j
)
+
∑
(
i
,
j
)
∈
A
p
i
j
d
i
j
}
s
.
t
.
{
t
i
j
≥
w
i
j
+
d
i
j
,
(
i
,
j
)
∈
A
x
j
−
x
i
≥
t
i
j
−
d
i
j
,
(
i
,
j
)
∈
A
d
i
j
≥
0
,
(
i
,
j
)
∈
A
\min\{15(x_8-x_1)+\sum\limits_{(i,j)\in A}q_{ij}(t_{ij}-d_{ij})+\sum\limits_{(i,j)\in A}p_{ij}d_{ij}\}\\ \begin{equation} s.t. \begin{cases} t_{ij}\geq w_{ij}+d_{ij},&(i,j)\in A\\ x_j-x_i\geq t_{ij}-d_{ij},&(i,j)\in A\\ d_{ij}\geq 0,&(i,j)\in A \end{cases} \end{equation}
min{15(x8−x1)+(i,j)∈A∑qij(tij−dij)+(i,j)∈A∑pijdij}s.t.⎩
⎨
⎧tij≥wij+dij,xj−xi≥tij−dij,dij≥0,(i,j)∈A(i,j)∈A(i,j)∈A
from scipy.optimize import linprog
import numpy as np
if __name__=="__main__":
f = np.zeros(18) # 由于需要求的最小表达式中有常数项，这里忽略了，如果要计算最终的费用需要加上
f[0] = -15
f[7] = 15
f[8:] = np.array([120-100,280-200,110-80,0-0,180-150,375-250,170-120,100-100,200-180,220-130])
A = np.array([[1,-1,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0],
[0,1,-1,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0],
[0,1,0,-1,0,0,0,0,0,0,-1,0,0,0,0,0,0,0],
[0,0,1,-1,0,0,0,0,0,0,0,-1,0,0,0,0,0,0],
[0,0,1,0,-1,0,0,0,0,0,0,0,-1,0,0,0,0,0],
[0,0,0,1,0,-1,0,0,0,0,0,0,0,-1,0,0,0,0],
[0,0,0,1,0,0,-1,0,0,0,0,0,0,0,-1,0,0,0],
[0,0,0,0,1,0,0,-1,0,0,0,0,0,0,0,-1,0,0],
[0,0,0,0,0,1,0,-1,0,0,0,0,0,0,0,0,-1,0],
[0,0,0,0,0,0,1,-1,0,0,0,0,0,0,0,0,0,-1]])
b = np.array([[-6,-9,-3,0,-7,-8,-2,-1,-4,-5]]).T
lb = np.zeros(18)
ub = np.array([None]*18)
ub[8:] = np.array([6-4,9-5,3-2,0-0,7-5,8-3,2-1,1-1,4-3,5-2])
bound = np.vstack((lb,ub)).T
res = linprog(f,A,b,bounds=bound)
print("需要花费的总天数为",res.x[7])
# output:需要花费的总天数为 27.0
扫一扫
8404
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
