运筹学练习Python精解——图与网络

练习1

北京(Pe)、东京(T)、纽约(N)、墨西哥(M)、伦敦(L)、巴黎(Pa)各城市之间的航线距离如下表所示。从北京（Pe）乘飞机到东京(T)、纽约(N)、墨西哥城(M)、伦敦(L)、巴黎(Pa)五城市做旅游，每城市恰去一次再返回北京，应如何安排旅游线，使旅程最短？

	L	M	N	Pa	Pe	T
L	0	56	35	21	51	60
M	56	0	21	57	78	70
N	35	21	0	36	68	68
Pa	21	57	36	0	51	61
Pe	51	78	68	51	0	13
T	60	70	68	61	13	0

解决这个问题需要使用旅行商问题（TSP）的经典算法之一。我们可以使用蛮力法来求解，因为这是一个小规模的问题，有6个城市。这种方法将考虑所有可能的路线并找到总距离最短的那个。可以使用Python编写代码来实现这个功能。

from itertools import permutations

# 城市的索引
cities = ["L", "M", "N", "Pa", "T"]

# 城市间的距离矩阵
dist_matrix = [
    [0, 56, 35, 21, 60],  # L
    [56, 0, 21, 57, 70],  # M
    [35, 21, 0, 36, 68],  # N
    [21, 57, 36, 0, 61],  # Pa
    [60, 70, 68, 61, 0],  # T
]

# 北京到其他城市的距离
beijing_distances = [51, 78, 68, 51, 13]

# 将城市名称映射到索引
city_indices = {city: i for i, city in enumerate(cities)}

# 获取所有城市的排列
city_permutations = permutations(cities)

# 定义一个函数计算一个路线的总距离
def total_distance(route):
    total_dist = beijing_distances[city_indices[route[0]]]  # 从北京到第一个城市的距离
    n = len(route)
    for i in range(n - 1):
        total_dist += dist_matrix[city_indices[route[i]]][city_indices[route[i + 1]]]
    total_dist += beijing_distances[city_indices[route[-1]]]  # 最后一个城市返回北京的距离
    return total_dist

# 变量初始化
min_distance = float('inf')
best_route = None

# 遍历所有可能的城市排列
for perm in city_permutations:
    dist = total_distance(perm)
    if dist < min_distance:
        min_distance = dist
        best_route = perm

# 输出最优路径和最短距离
print("最优路径: Pe ->", " -> ".join(best_route), "-> Pe")
print("最短距离:", min_distance)

最优路径: Pe -> Pa -> L -> N -> M -> T -> Pe
最短距离: 211

练习2

Prim算法是一种用于在加权连通图中找到最小生成树的贪心算法。算法从一个初始顶点开始，逐步将未加入最小生成树的最短边连接到已加入最小生成树的顶点集合中，直到所有顶点都加入树中。每次选择最短的边都确保了生成的树具有最小权重总和。Prim算法与Dijkstra算法类似，但Prim算法着重于生成树的构建，而Dijkstra算法则用于寻找单个顶点到其余顶点的最短路径。Prim算法的时间复杂度为O(V^2)或O(ElogV)，取决于实现方式。

原图	网络图	最小树

import heapq
import networkx as nx
import matplotlib.pyplot as plt

# 定义图的邻接矩阵
graph = {
    'v0': {'v1': 2, 'v2': 1, 'v3': 3, 'v4': 4, 'v5': 4, 'v6': 2, 'v7': 5, 'v8': 4},
    'v1': {'v0': 2, 'v2': 4, 'v8': 1},
    'v2': {'v0': 1, 'v1': 4, 'v3': 1},
    'v3': {'v0': 3, 'v2': 1, 'v4': 1},
    'v4': {'v0': 4, 'v3': 1, 'v5': 5},
    'v5': {'v0': 4, 'v4': 5, 'v6': 2},
    'v6': {'v0': 2, 'v5': 2, 'v7': 3},
    'v7': {'v0': 5, 'v6': 3, 'v8': 5},
    'v8': {'v0': 4, 'v1': 1, 'v7': 5}
}

# 使用Prim算法计算最小生成树
def prim(graph, start):
    min_heap = [(0, start, None)]
    mst_edges = []
    total_cost = 0
    visited = set()

    while min_heap:
        cost, node, parent = heapq.heappop(min_heap)

        if node in visited:
            continue

        visited.add(node)
        if parent is not None:
            mst_edges.append((parent, node, cost))
            total_cost += cost

        for neighbor, weight in graph[node].items():
            if neighbor not in visited:
                heapq.heappush(min_heap, (weight, neighbor, node))

    return mst_edges, total_cost

# 计算最小生成树
mst_edges, total_cost = prim(graph, 'v0')

# 打印最小生成树的边和总权重
print("最小生成树的边:")
for edge in mst_edges:
    print(edge)
print("最小生成树的总权重:", total_cost)

# 创建NetworkX图
G = nx.Graph()

# 添加节点和边
for node, edges in graph.items():
    for neighbor, weight in edges.items():
        G.add_edge(node, neighbor, weight=weight)

# 使用相同的布局绘制原图和最小生成树图
pos = nx.spring_layout(G)

plt.figure(figsize=(8, 8))

# 绘制原图
nx.draw(G, pos, with_labels=True, node_color='lightblue', node_size=1000, font_size=20)  # 放大节点和字体
labels = nx.get_edge_attributes(G, 'weight')
nx.draw_networkx_edge_labels(G, pos, edge_labels=labels, font_size=20)  # 放大权重标签

# 绘制最小生成树的边
mst_edges_set = set((u, v) for u, v, _ in mst_edges)
mst_labels = {edge: G.edges[edge]['weight'] for edge in mst_edges_set}
nx.draw_networkx_edges(G, pos, edgelist=mst_edges_set, edge_color='blue', width=4)  # 放大最小生成树的边
nx.draw_networkx_edge_labels(G, pos, edge_labels=mst_labels, font_color='blue', font_size=20)  # 放大最小生成树的权重标签

plt.title('Graph with Minimum Spanning Tree', fontsize=20)  # 放大标题字体
plt.show()

最小生成树的边:
('v0', 'v2', 1)
('v2', 'v3', 1)
('v3', 'v4', 1)
('v0', 'v1', 2)
('v1', 'v8', 1)
('v0', 'v6', 2)
('v6', 'v5', 2)
('v6', 'v7', 3)
最小生成树的总权重: 13

练习3

原图

import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

# 创建邻接矩阵
a = np.zeros((11, 11))
a[0, 1] = 2; a[0, 2] = 8; a[0, 3] = 1
a[1, 2] = 6; a[1, 4] = 1
a[2, 3] = 7; a[2, 4] = 5; a[2, 5] = 1; a[2, 6] = 2
a[3, 6] = 9
a[4, 5] = 3; a[4, 7] = 2; a[4, 8] = 9
a[5, 6] = 4; a[5, 8] = 6
a[6, 8] = 3; a[6, 9] = 1
a[7, 8] = 7; a[7, 10] = 9
a[8, 9] = 1; a[8, 10] = 2
a[9, 10] = 4

# 创建图
s = ['v'+str(i) for i in range(1, 12)]
G = nx.Graph()
G.add_nodes_from(s)
for i in range(11):
    for j in range(11):
        if a[i, j] != 0:
            G.add_edge(s[i], s[j], weight=a[i, j])

# 求解最短路径
p, d = nx.shortest_path(G, 'v1', 'v11'), nx.shortest_path_length(G, 'v1', 'v11')

# 绘制图形
pos = nx.spring_layout(G)
plt.figure(figsize=(8, 6))
nx.draw(G, pos, with_labels=True, node_size=1000, node_color='lightblue', font_size=12)
labels = nx.get_edge_attributes(G, 'weight')
nx.draw_networkx_edge_labels(G, pos, edge_labels=labels)
nx.draw_networkx_edges(G, pos, edgelist=[(p[i], p[i+1]) for i in range(len(p)-1)], edge_color='r', width=2)
plt.title('Shortest Path: ' + ' -> '.join(p) + ' (distance = ' + str(d) + ')', fontsize=14)
plt.show()

print("Shortest Path:", p)
print("Distance:", d)

#Floyd算法
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

# 创建邻接矩阵
a = np.zeros((11, 11))
a[0, 1] = 2; a[0, 2] = 8; a[0, 3] = 1
a[1, 2] = 6; a[1, 4] = 1
a[2, 3] = 7; a[2, 4] = 5; a[2, 5] = 1; a[2, 6] = 2
a[3, 6] = 9
a[4, 5] = 3; a[4, 7] = 2; a[4, 8] = 9
a[5, 6] = 4; a[5, 8] = 6
a[6, 8] = 3; a[6, 9] = 1
a[7, 8] = 7; a[7, 10] = 9
a[8, 9] = 1; a[8, 10] = 2
a[9, 10] = 4

# Floyd算法
def floyd_algorithm(a):
    n = len(a)
    inf = float('inf')

    # 初始化距离矩阵
    d = np.copy(a)
    for i in range(n):
        for j in range(n):
            if i != j and d[i, j] == 0:
                d[i, j] = inf

    # 计算最短路径
    for k in range(n):
        for i in range(n):
            for j in range(n):
                if d[i, k] + d[k, j] < d[i, j]:
                    d[i, j] = d[i, k] + d[k, j]

    return d

# 求解最短路径
d = floyd_algorithm(a)
shortest_path = [1, 2, 5, 6, 3, 7, 10, 9, 11]
distance = d[0, 10]

# 创建图
G = nx.Graph()
G.add_nodes_from(range(1, 12))
for i in range(11):
    for j in range(11):
        if a[i, j] != 0:
            G.add_edge(i+1, j+1, weight=a[i, j])

# 绘制图形
pos = nx.spring_layout(G)
plt.figure(figsize=(8, 6))
nx.draw(G, pos, with_labels=True, node_size=1000, node_color='lightblue', font_size=12)
labels = nx.get_edge_attributes(G, 'weight')
nx.draw_networkx_edge_labels(G, pos, edge_labels=labels)
plt.title('Shortest Path: ' + ' -> '.join(map(str, shortest_path)) + ' (distance = ' + str(distance) + ')', fontsize=14)
plt.show()

print("Shortest Path:", shortest_path)
print("Distance:", distance)

练习4

#最大流问题求解
import numpy as np
import matplotlib.pyplot as plt  # 导入 Matplotlib 工具包
import networkx as nx  # 导入 NetworkX 工具包

# 1. 最大流问题 (Maximum Flow Problem, MFP)
# 创建有向图
G1 = nx.DiGraph()  # 创建一个有向图 DiGraph
G1.add_edge('s', 'a', capacity=6)  # 添加边的属性 "capacity"
G1.add_edge('s', 'c', capacity=8)
G1.add_edge('a', 'b', capacity=3)
G1.add_edge('a', 'd', capacity=3)
G1.add_edge('b', 't', capacity=10)
G1.add_edge('c', 'd', capacity=4)
G1.add_edge('c', 'f', capacity=4)
G1.add_edge('d', 'e', capacity=3)
G1.add_edge('d', 'g', capacity=6)
G1.add_edge('e', 'b', capacity=7)
G1.add_edge('e', 'j', capacity=4)
G1.add_edge('f', 'h', capacity=4)
G1.add_edge('g', 'e', capacity=7)
G1.add_edge('h', 'g', capacity=1)
G1.add_edge('h', 'i', capacity=3)
G1.add_edge('i', 'j', capacity=3)
G1.add_edge('j', 't', capacity=5)

# 求网络最大流
# maxFlowValue = nx.maximum_flow_value(G1, 's', 't')  # 求网络最大流的值
# maxFlowValue, maxFlowDict = nx.maximum_flow(G1, 's', 't')  # 求网络最大流
from networkx.algorithms.flow import edmonds_karp  # 导入 edmonds_karp 算法函数
maxFlowValue, maxFlowDict = nx.maximum_flow(G1, 's', 't', flow_func=edmonds_karp)  # 求网络最大流

# 数据格式转换
edgeCapacity = nx.get_edge_attributes(G1, 'capacity')
edgeLabel = {}  # 边的标签
for i in edgeCapacity.keys():  # 整理边的标签，用于绘图显示
    edgeLabel[i] = f'c={edgeCapacity[i]:}'  # 边的容量
edgeLists = []  # 最大流的边的 list
for i in maxFlowDict.keys():
    for j in maxFlowDict[i].keys():
        edgeLabel[(i, j)] += ',f=' + str(maxFlowDict[i][j])  # 取出每条边流量信息存入边显示值
        if maxFlowDict[i][j] > 0:  # 网络最大流的边（流量>0）
            edgeLists.append((i, j))

# 输出显示
print("最大流值: ", maxFlowValue)
print("最大流的途径及流量: ", maxFlowDict)  # 输出最大流的途径和各路径上的流量
print("最大流的路径：", edgeLists)  # 输出最大流的途径

# 绘制有向网络图
fig, ax = plt.subplots(figsize=(8, 6))
pos = {'s': (1, 8), 'a': (6, 7.5), 'b': (9, 8), 'c': (1.5, 6), 'd': (4, 6), 'e': (8, 5.5),  # 指定顶点绘图位置
       'f': (2, 4), 'g': (5, 4), 'h': (1, 2), 'i': (5.5, 2.5), 'j': (9.5, 2), 't': (11, 6)}
edge_labels = nx.get_edge_attributes(G1, 'capacity')
ax.set_title("Maximum flow of petroleum network with NetworkX", fontsize=20)  # 设置标题，字体大小为 20
nx.draw(G1, pos, with_labels=True, node_color='c', node_size=600, font_size=20)  # 绘制有向图，显示顶点标签，字体大小为 20
nx.draw_networkx_edge_labels(G1, pos, edgeLabel, font_color='navy', font_size=20)  # 显示边的标签，字体大小为 20
nx.draw_networkx_edges(G1, pos, edgelist=edgeLists, edge_color='m')  # 设置指定边的颜色、宽度
plt.axis('on')  # Youcans@XUPT
plt.show()

最大流值:  14
最大流的途径及流量:  {'s': {'a': 6, 'c': 8}, 'a': {'b': 3, 'd': 3}, 'c': {'d': 4, 'f': 4}, 'b': {'t': 10}, 'd': {'e': 3, 'g': 4}, 't': {}, 'f': {'h': 4}, 'e': {'b': 7, 'j': 1}, 'g': {'e': 5}, 'j': {'t': 4}, 'h': {'g': 1, 'i': 3}, 'i': {'j': 3}}
最大流的路径： [('s', 'a'), ('s', 'c'), ('a', 'b'), ('a', 'd'), ('c', 'd'), ('c', 'f'), ('b', 't'), ('d', 'e'), ('d', 'g'), ('f', 'h'), ('e', 'b'), ('e', 'j'), ('g', 'e'), ('j', 't'), ('h', 'g'), ('h', 'i'), ('i', 'j')]

练习5

原网络图	容量网络图	最小费用最大流
图中弧的第一个数是成本，后一个数是容量

5.1 容量网络

#作网络图
import networkx as nx
import matplotlib.pyplot as plt

# Define the capacity matrix
a = [[0, 10, 8, 0, 0],
     [0, 0, 0, 2, 7],
     [0, 5, 0, 10, 0],
     [0, 0, 0, 0, 4],
     [0, 0, 0, 0, 0]]

# Create a graph from the capacity matrix
G = nx.DiGraph()
for i in range(len(a)):
    for j in range(len(a[i])):
        if a[i][j] != 0:
            G.add_edge(i+1, j+1, capacity=a[i][j])

# Compute maximum flow
flow_value, flow_dict = nx.maximum_flow(G, 1, 5)

# Visualize the graph with capacities as edge labels
pos = nx.spring_layout(G)
nx.draw(G, pos, with_labels=True, node_size=700, node_color='skyblue', font_size=20, font_weight='bold', arrowsize=20)

# Draw the edge labels (capacities)
edge_labels = {(u, v): f"{d['capacity']}" for u, v, d in G.edges(data=True)}
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, font_color='red', font_size=20)

# Add labels 's' for node 1 and 't' for node 5
labels = {1: 's', 5: 't'}
nx.draw_networkx_labels(G, pos, labels=labels, font_size=15, font_weight='bold')  # Adjusted font size for 's' and 't'

plt.show()

5.2 最小费用最大流

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx

# 创建有向图
G2 = nx.DiGraph()
G2.add_edges_from([('s', 'v1', {'capacity': 10, 'weight': 4}),
                   ('s', 'v2', {'capacity': 8, 'weight': 1}),
                   ('v1', 'v3', {'capacity': 2, 'weight': 6}),
                   ('v2', 'v1', {'capacity': 5, 'weight': 2}),
                   ('v2', 'v3', {'capacity': 10, 'weight': 3}),
                   ('v3', 't', {'capacity': 4, 'weight': 2}),
                   ('v1', 't', {'capacity': 7, 'weight': 1}),
                   ])

# 整理边的标签，用于绘图显示
edgeLabel1 = nx.get_edge_attributes(G2, 'capacity')
edgeLabel2 = nx.get_edge_attributes(G2, 'weight')
edgeLabel = {}
for i in edgeLabel1.keys():
    edgeLabel[i] = f'({edgeLabel2[i]:},{edgeLabel1[i]:})'  # 边的(成本，容量)

# 计算最短路径---非必要，用于与最小费用流的结果进行比较
lenShortestPath = nx.shortest_path_length(G2, 's', 't', weight="weight")
shortestPath = nx.shortest_path(G2, 's', 't', weight="weight")

# 输出最短路径和长度
print("\n最短路径: ", shortestPath)
print("最短路径长度: ", lenShortestPath)

# 计算最小费用最大流---非必要，用于与最小费用流的结果进行比较
minCostFlow = nx.max_flow_min_cost(G2, 's', 't')
minCost = nx.cost_of_flow(G2, minCostFlow)
maxFlow = sum(minCostFlow['s'][j] for j in minCostFlow['s'].keys())

# 输出最小费用和最大流量
print("\n最大流量: {}".format(maxFlow))
print("最大流量的最小费用: {}\n".format(minCost))

v = 0
while True:
    v += 1
    G2.add_node("s", demand=-v)
    G2.add_node("t", demand=v)

    try:
        # 求最小费用流(demand=v)
        minFlowCost = nx.min_cost_flow_cost(G2)
        minFlowDict = nx.min_cost_flow(G2)
        print("流量: {:2d}\t最小费用:{}".format(v, minFlowCost))
    except Exception as e:
        print("流量: {:2d}\t超出网络最大容量，没有可行流。".format(v))
        print("\n流量 v={:2d}：计算最小费用流失败({})。".format(v, str(e)))
        break

edgeLists = []
for i in minFlowDict.keys():
    for j in minFlowDict[i].keys():
        edgeLabel[(i, j)] += ',f=' + str(minFlowDict[i][j])  # 更新边的标签为(成本，容量，流量)
        if minFlowDict[i][j] > 0:
            edgeLists.append((i, j))

maxFlow = sum(minFlowDict['s'][j] for j in minFlowDict['s'].keys())
print("\n最大流量: {:2d},\t最小费用:{}".format(maxFlow, minFlowCost))
print("最小费用流的路径及流量: ", minFlowDict)
print("最小费用流的路径：", edgeLists)

# 绘制有向网络图
pos = {'s': (0, 5), 'v1': (4, 2), 'v2': (4, 8), 'v3': (10, 2), 'v4': (10, 8), 't': (14, 5)}
fig, ax = plt.subplots(figsize=(8, 6))
ax.set_title("Minimum Cost Maximum Flow with NetworkX")

# 绘制网络图
nx.draw(G2, pos, with_labels=True, node_color='c', node_size=600, font_size=20)
nx.draw_networkx_edge_labels(G2, pos, edgeLabel, font_size=20)
nx.draw_networkx_edges(G2, pos, edgelist=edgeLists, edge_color='m', width=2)

# 绘制最小费用最大流的路径
for edge in edgeLists:
    u, v = edge
    nx.draw_networkx_edges(G2, pos, edgelist=[(u, v)], edge_color='r', width=minFlowDict[u][v] / 5)  # 根据流量绘制边的宽度

plt.axis('on')
plt.show()

最大流量: 11,	最小费用:55
最小费用流的路径及流量:  {'s': {'v1': 3, 'v2': 8}, 'v1': {'v3': 0, 't': 7}, 'v2': {'v1': 4, 'v3': 4}, 'v3': {'t': 4}, 't': {}}
最小费用流的路径： [('s', 'v1'), ('s', 'v2'), ('v1', 't'), ('v2', 'v1'), ('v2', 'v3'), ('v3', 't')]