尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

用Python复现数学建模经典赛题:Floyd算法与0-1规划实战

用Python复现数学建模经典赛题:Floyd算法与0-1规划实战 1. 项目背景与核心挑战十多年前2011年的全国大学生数学建模竞赛B题《交巡警服务平台的设置与调度》可以说是一道经典的优化问题也是很多同学接触运筹学、图论和编程求解的“启蒙题”。这道题的核心是要求参赛者在一个城市的交通网络图上为交巡警服务平台规划合理的管辖范围并设计一个在发生突发事件时从20个平台中快速调度13辆警车前往13个出入城区的路口进行全城封锁的方案。题目给出了城市92个节点的坐标、相邻节点间的道路长度以及20个平台的初始位置。当年很多优秀论文的求解思路都离不开几个关键工具用Floyd算法计算任意两点间的最短路径用0-1整数规划模型来精确描述平台管辖分配和警车调度问题然后用Lingo这类专业的优化软件来求解模型。对于当时的大学生来说能熟练运用这些工具并写出逻辑清晰的论文已经非常了不起了。然而时过境迁技术栈发生了翻天覆地的变化。今天我们有了更强大、更易用的工具。Python凭借其丰富的数据科学生态如NumPy, Pandas, SciPy和优化库如PuLP, OR-Tools成为了解决此类问题的首选。Excel则从单纯的数据展示工具进化成了与Python无缝衔接的数据预处理和结果可视化的得力助手。复现这道题不再是为了追求一个标准答案而是为了用现代的、可复现的技术流水线去重新审视和解决一个经典的运筹优化问题。这个过程本身就是对算法思想、建模技巧和工程实现能力的一次综合演练。2. 数据准备与预处理从原始数据到可计算网络任何建模工作的第一步都是处理数据。原题数据通常以文本或表格形式给出节点坐标和邻接关系。我们的目标是将这些原始数据转化为Python中易于操作的结构并为后续的图论计算做好准备。2.1 原始数据解析与清洗题目数据通常包含三部分节点信息92个节点的编号、横纵坐标单位毫米基于地图比例尺。边信息节点之间的连接关系即道路以及对应的实际长度。平台信息20个服务平台的初始设置节点编号。首先我们需要将这些数据数字化。虽然当年可能用Matlab或手工处理但现在用Pandas读取是更高效的方式。假设我们有一个data.xlsx的Excel文件里面按照上述结构整理了三个工作表。import pandas as pd # 读取节点数据 nodes_df pd.read_excel(data.xlsx, sheet_name节点) # 假设列名为[节点编号, 横坐标, 纵坐标] print(nodes_df.head()) # 读取边数据 edges_df pd.read_excel(data.xlsx, sheet_name道路) # 假设列名为[起点, 终点, 长度] print(edges_df.head()) # 读取平台数据 platforms_df pd.read_excel(data.xlsx, sheet_name平台) # 假设列名为[平台编号, 所在节点] print(platforms_df.head())数据清洗可能包括检查是否有重复的边如既记录了A-B又记录了B-A、是否存在孤立的节点、长度是否为非负数等。对于本问题由于是无向图我们需要确保每条道路在计算时被正确视为双向通行。# 简单的数据清洗示例 # 确保长度为正 if (edges_df[长度] 0).any(): print(警告存在非正长度的道路请检查数据。) # 通常可以过滤掉这些异常边 edges_df edges_df[edges_df[长度] 0] # 由于是无向图我们可以复制一份反向的边确保邻接矩阵的对称性但Floyd算法通常用有向图矩阵表示无向图对称初始化即可2.2 构建图论模型的基础数据结构接下来我们需要将DataFrame形式的数据转化为图论算法如Floyd算法可以直接使用的数据结构。最常用的就是邻接矩阵。对于一个有n个节点本题n92的图我们初始化一个n x n的矩阵矩阵中的值dist[i][j]表示从节点i到节点j的最短距离。初始时如果i和j之间有直接相连的边则dist[i][j]等于该边长度否则dist[i][j]被设为一个非常大的数代表无穷大即不可直接到达对角线元素dist[i][i]设为0。这里有一个关键细节节点编号通常是从1开始的而Python的列表索引是从0开始的。我们需要建立一个映射关系。import numpy as np n_nodes 92 # 创建一个节点编号到矩阵索引0-91的映射字典 node_id_to_index {node_id: idx for idx, node_id in enumerate(nodes_df[节点编号])} index_to_node_id {idx: node_id for node_id, idx in node_id_to_index.items()} # 初始化邻接矩阵用一个大数如np.inf表示无穷远 dist_matrix np.full((n_nodes, n_nodes), np.inf) # 将对角线设为0 np.fill_diagonal(dist_matrix, 0) # 填充直接相连的边 for _, row in edges_df.iterrows(): u node_id_to_index[row[起点]] v node_id_to_index[row[终点]] length row[长度] # 无向图双向距离相同 dist_matrix[u][v] length dist_matrix[v][u] length print(f邻接矩阵形状{dist_matrix.shape}) print(f示例节点1到节点2的距离如果相连: {dist_matrix[node_id_to_index[1], node_id_to_index[2]]})至此我们得到了一个包含92个节点、基于实际道路长度的加权邻接矩阵这是所有后续优化计算的基础。3. 核心算法一Floyd算法求解全局最短路径在分配平台管辖范围和调度警车时核心依据是“距离”或“时间”这需要我们知道图中任意两个节点之间的最短路径长度。Floyd-Warshall算法正是解决全源最短路径问题的经典动态规划算法。3.1 Floyd算法原理与实现Floyd算法的思想非常巧妙它通过考虑逐步引入中间节点k来更新任意两点i和j之间的最短距离。状态转移方程为dist[i][j] min(dist[i][j], dist[i][k] dist[k][j])其中k是中间节点。我们需要对所有的k从0到n-1进行循环。def floyd_warshall(dist): 使用Floyd-Warshall算法计算所有节点对之间的最短距离。 参数 dist: 初始化的邻接矩阵dist[i][j]表示i到j的直接距离无边则为inf。 返回: 更新后的最短距离矩阵。 n dist.shape[0] # 复制一份避免修改原矩阵 d dist.copy() # 三重循环是算法的核心 for k in range(n): # 中间节点 for i in range(n): # 起点 # 一个小优化如果i到k的距离是无穷大则跳过j的循环 if d[i, k] np.inf: continue for j in range(n): # 终点 # 如果i-k和k-j都是可达的则尝试更新i-j的距离 if d[k, j] np.inf: continue new_dist d[i, k] d[k, j] if new_dist d[i, j]: d[i, j] new_dist return d # 计算全局最短路径矩阵 shortest_path_matrix floyd_warshall(dist_matrix) print(全局最短路径矩阵计算完成。) print(f示例节点1到节点92的最短距离: {shortest_path_matrix[node_id_to_index[1], node_id_to_index[92]]})注意原始的Floyd算法时间复杂度是O(n³)对于n92来说完全在可接受范围内约92^3 ≈ 78万次操作。但如果节点数上万就需要考虑更高效的算法如多次Dijkstra算法或使用优化库了。这里的实现是一个清晰的演示在实际复现中也可以直接使用scipy.sparse.csgraph.floyd_warshall函数它经过优化可能更快。3.2 最短路径矩阵的应用与验证得到shortest_path_matrix后我们就拥有了解决后续所有优化问题的“距离字典”。例如查询平台节点p到任意案件发生节点a的最短距离shortest_path_matrix[index_p, index_a]为每个节点分配管辖平台时就是找距离它最短的那个平台。我们可以快速验证一下算法的正确性比如检查从某个平台到所有节点的距离是否都是有限值即全图连通以及对称性无向图距离矩阵应对称。# 检查连通性是否存在无穷大值除了对角线 non_diag shortest_path_matrix[~np.eye(shortest_path_matrix.shape[0], dtypebool)] if np.isinf(non_diag).any(): print(警告图中存在不连通的节点对) else: print(图是连通的所有节点对之间都有路径。) # 检查对称性由于计算误差用近似相等判断 if not np.allclose(shortest_path_matrix, shortest_path_matrix.T): print(警告最短路径矩阵不对称可能计算有误或原图是有向图。) else: print(最短路径矩阵对称性良好。)这个最短路径矩阵是整个项目的基石后面的0-1规划模型将重度依赖它提供的距离数据。4. 核心算法二0-1整数规划建模与求解有了距离数据我们就可以用数学规划的语言来精确描述题目中的两个优化问题平台管辖分配和警车调度。这两个问题本质上都是组合优化问题非常适合用0-1整数规划决策变量只能取0或1来建模。4.1 问题一服务平台管辖范围分配建模目标为20个服务平台分配管辖节点使得“所有节点到其管辖平台的最大距离”最小化。这是一个典型的最小化最大距离问题也称为“中心点”或“最小覆盖半径”问题。决策变量x_{ij}: 0-1变量。如果节点i共92个由平台j共20个管辖则为1否则为0。R: 一个连续变量表示所有分配方案中的最大距离即覆盖半径。约束条件每个节点必须且只能被一个平台管辖对于每一个节点i∑_j x_{ij} 1。一个平台可以管辖多个节点但没有上限约束原题如此。距离约束对于每一对(i, j)如果x_{ij}1那么节点i到平台j的距离d_{ij}必须小于等于R。这是一个逻辑约束需要线性化。标准方法是d_{ij} * x_{ij} R。因为当x_{ij}0时不等式自动成立当x_{ij}1时要求d_{ij} R。变量定义域x_{ij} ∈ {0, 1},R 0。目标函数最小化R。这个模型可以直接用Python的优化库来构建和求解。我们以流行的PuLP库为例。import pulp # 准备数据 n_nodes 92 platform_indices [node_id_to_index[pid] for pid in platforms_df[所在节点]] # 20个平台的索引 # 距离字典方便查询 dist_dict {(i, j): shortest_path_matrix[i, j] for i in range(n_nodes) for j in platform_indices} # 创建问题 prob1 pulp.LpProblem(Platform_Allocation_MinMax, pulp.LpMinimize) # 创建决策变量 x pulp.LpVariable.dicts(x, ((i, j) for i in range(n_nodes) for j in platform_indices), lowBound0, upBound1, catBinary) R pulp.LpVariable(R, lowBound0, catContinuous) # 设置目标函数 prob1 R # 添加约束 # 1. 每个节点只能由一个平台管辖 for i in range(n_nodes): prob1 pulp.lpSum(x[i, j] for j in platform_indices) 1 # 2. 距离约束线性化后的形式 for i in range(n_nodes): for j in platform_indices: prob1 dist_dict[(i, j)] * x[i, j] R # 求解问题 solver pulp.PULP_CBC_CMD(msgFalse) # 使用CBC求解器不输出求解日志 prob1.solve(solver) # 输出结果 print(f问题一状态: {pulp.LpStatus[prob1.status]}) print(f最小化最大距离 R {pulp.value(R):.2f}) # 解析分配方案 allocation {} for i in range(n_nodes): for j in platform_indices: if pulp.value(x[i, j]) 1: allocation[index_to_node_id[i]] index_to_node_id[j] break print(f分配方案已计算例如节点1由平台{allocation[1]}管辖。)4.2 问题二警车调度优化建模目标从20个平台中选出13个各派出一辆警车前往13个指定的出入城路口进行封锁使得“警车行驶总距离最短”。这是一个指派问题或二分图最小权匹配问题。决策变量y_{kj}: 0-1变量。如果第k个封锁任务对应第k个路口共13个由平台j共20个的警车执行则为1否则为0。约束条件每个封锁路口必须有一辆警车前往对于每一个路口k∑_j y_{kj} 1。每个平台最多派出一辆警车因为只有13个任务从20个中选13个对于每一个平台j∑_k y_{kj} 1。正好有13辆警车被派出∑_k ∑_j y_{kj} 13。目标函数最小化总行驶距离 ∑_k ∑_j (d_{kj} * y_{kj})其中d_{kj}是平台j到路口k的最短距离。# 假设我们有一个列表包含了13个需要封锁的路口节点编号 blockade_nodes [假设的13个路口节点ID列表] blockade_indices [node_id_to_index[nid] for nid in blockade_nodes] # 创建问题 prob2 pulp.LpProblem(Police_Dispatch_MinSum, pulp.LpMinimize) # 创建决策变量 y pulp.LpVariable.dicts(y, ((k_idx, j_idx) for k_idx in blockade_indices for j_idx in platform_indices), catBinary) # 设置目标函数最小化总距离 prob2 pulp.lpSum(shortest_path_matrix[k_idx, j_idx] * y[k_idx, j_idx] for k_idx in blockade_indices for j_idx in platform_indices) # 添加约束 # 1. 每个路口必须有一辆警车 for k_idx in blockade_indices: prob2 pulp.lpSum(y[k_idx, j_idx] for j_idx in platform_indices) 1 # 2. 每个平台最多派出一辆警车 for j_idx in platform_indices: prob2 pulp.lpSum(y[k_idx, j_idx] for k_idx in blockade_indices) 1 # 3. 正好派出13辆警车此约束可由前两个约束推导但显式写出更清晰 prob2 pulp.lpSum(y[k_idx, j_idx] for k_idx in blockade_indices for j_idx in platform_indices) 13 # 求解 prob2.solve(solver) print(f\n问题二状态: {pulp.LpStatus[prob2.status]}) print(f最小化总距离 {pulp.value(prob2.objective):.2f}) # 解析调度方案 dispatch_plan [] for k_idx in blockade_indices: for j_idx in platform_indices: if pulp.value(y[k_idx, j_idx]) 1: dispatch_plan.append((index_to_node_id[k_idx], index_to_node_id[j_idx], shortest_path_matrix[k_idx, j_idx])) break print(调度方案路口-平台距离) for plan in dispatch_plan: print(f 路口{plan[0]} - 平台{plan[1]}, 距离{plan[2]:.2f})通过以上两个模型我们就能得到理论上最优的分配和调度方案。对比当年优秀论文的结果可以用现代工具验证其正确性或发现更优解。5. 工程实现与可视化从结果到洞察得到数学上的最优解只是第一步。作为一个完整的复现项目我们需要将结果清晰、直观地呈现出来并确保整个流程的工程化使其可重复、可验证。5.1 结果分析与Excel输出将Python的计算结果导出到Excel可以方便地进行进一步分析、制作报告或与他人协作。Pandas的to_excel功能非常强大。# 1. 将问题一的分配方案转为DataFrame allocation_list [] for node_id, plat_id in allocation.items(): distance shortest_path_matrix[node_id_to_index[node_id], node_id_to_index[plat_id]] allocation_list.append({节点编号: node_id, 管辖平台: plat_id, 最短距离: distance}) allocation_df pd.DataFrame(allocation_list) # 计算每个平台的管辖节点数和工作量总距离或最大距离 platform_summary allocation_df.groupby(管辖平台).agg( 管辖节点数(节点编号, count), 最远距离(最短距离, max), 平均距离(最短距离, mean) ).reset_index() # 2. 将问题二的调度方案转为DataFrame dispatch_df pd.DataFrame(dispatch_plan, columns[封锁路口, 出警平台, 行驶距离]) # 3. 写入Excel的不同工作表 with pd.ExcelWriter(2011B题_求解结果.xlsx, engineopenpyxl) as writer: allocation_df.to_excel(writer, sheet_name平台管辖分配明细, indexFalse) platform_summary.to_excel(writer, sheet_name平台工作量统计, indexFalse) dispatch_df.to_excel(writer, sheet_name警车调度方案, indexFalse) # 还可以把原始数据也放进去 nodes_df.to_excel(writer, sheet_name原始节点数据, indexFalse) edges_df.to_excel(writer, sheet_name原始道路数据, indexFalse) print(结果已成功导出至 2011B题_求解结果.xlsx)在Excel中我们可以利用条件格式、图表等功能直观展示哪个平台管辖范围最大、哪个路口距离最远、调度方案是否均衡等。5.2 网络图可视化一图胜千言。使用networkx和matplotlib或更现代的plotly可以将网络、平台、分配方案和调度路线可视化。import networkx as nx import matplotlib.pyplot as plt # 创建图对象 G nx.Graph() # 添加节点附带位置属性 for _, row in nodes_df.iterrows(): node_id row[节点编号] G.add_node(node_id, pos(row[横坐标], row[纵坐标])) # 添加边 for _, row in edges_df.iterrows(): G.add_edge(row[起点], row[终点], weightrow[长度]) # 获取位置字典 pos nx.get_node_attributes(G, pos) plt.figure(figsize(15, 10)) # 1. 绘制所有节点和边浅灰色背景网络 nx.draw_networkx_nodes(G, pos, node_size20, node_colorlightgray, alpha0.6) nx.draw_networkx_edges(G, pos, edge_colorgray, alpha0.4, width1) # 2. 高亮平台节点 platform_nodes platforms_df[所在节点].tolist() nx.draw_networkx_nodes(G, pos, nodelistplatform_nodes, node_size200, node_colorred, label服务平台) # 3. 高亮封锁路口节点 blockade_nodes_list [plan[0] for plan in dispatch_plan] # 从调度方案中提取 nx.draw_networkx_nodes(G, pos, nodelistblockade_nodes_list, node_size150, node_colorblue, label封锁路口) # 4. 用不同颜色标记不同平台的管辖范围示例为前5个平台着色 colors [green, orange, purple, brown, pink] for idx, plat_id in enumerate(platform_nodes[:5]): # 找出由该平台管辖的所有节点 covered_nodes [node_id for node_id, p_id in allocation.items() if p_id plat_id] nx.draw_networkx_nodes(G, pos, nodelistcovered_nodes, node_size30, node_colorcolors[idx % len(colors)], alpha0.7) # 5. 绘制调度路线从平台到路口 for plan in dispatch_plan: # 注意这里绘制的是直接连线并非实际最短路径。若要绘制真实路径需用Floyd算法记录前驱节点并回溯。 nx.draw_networkx_edges(G, pos, edgelist[(plan[1], plan[0])], edge_colorred, width2, styledashed, alpha0.8) plt.title(2011国赛B题交巡警服务平台设置与调度方案可视化) plt.legend(scatterpoints1) plt.axis(off) # 关闭坐标轴 plt.tight_layout() plt.savefig(solution_visualization.png, dpi300) plt.show()可视化能立刻让我们看到平台分布是否均匀、管辖范围是否紧凑、调度路线是否交叉严重等问题这是纯数据表格无法提供的直观感受。5.3 项目代码的组织与复现性保障为了让这个复现项目真正具有价值而不仅仅是一次性脚本良好的代码组织至关重要。2011_MCM_B_Reproduction/ ├── data/ │ ├── raw/ # 存放原始题目数据文件 │ └── processed/ # 存放清洗后的数据文件 ├── src/ │ ├── 01_data_preprocessing.py │ ├── 02_floyd_algorithm.py │ ├── 03_model_building.py # 包含两个0-1规划模型 │ ├── 04_visualization.py │ └── utils.py # 公共函数如节点映射 ├── results/ │ ├── excel/ # 输出的Excel结果 │ └── figures/ # 生成的图表 ├── requirements.txt # 项目依赖库列表 ├── README.md # 项目说明包括问题重述、运行方法 └── main.py # 主运行脚本按顺序调用各模块在requirements.txt中写明依赖pandas1.3.0 numpy1.21.0 pulp2.6.0 networkx2.6.0 matplotlib3.4.0 openpyxl3.0.0在README.md中详细说明如何安装环境、运行代码、解读结果。这样任何人拿到这个项目都能一键复现全部过程和结果这正是现代数据科学项目的标准做法。6. 踩坑实录与进阶思考复现一个十多年前的赛题绝不仅仅是把答案算出来。在这个过程中你会遇到许多当年论文可能一笔带过但实际编程中却至关重要的问题。6.1 Floyd算法中的“无穷大”与数值稳定性在初始化邻接矩阵时我们用np.inf代表无穷大。但在Floyd算法的循环中需要进行加法比较d[i, k] d[k, j]。如果d[i, k]或d[k, j]是np.inf那么它们的和也是np.inf与另一个np.inf比较在Python中是安全的。但为了效率和避免不必要的计算我在代码中加入了if判断来跳过无效循环。这是一个常见的微优化。另一个坑是数值精度。道路长度可能是浮点数多次累加后可能产生微小的误差。在检查对称性时我使用了np.allclose而不是直接判断相等。在后续的优化模型中距离作为系数输入这些微小误差通常不影响0-1规划的最优解但严谨起见需要注意。6.2 0-1规划模型构建的陷阱问题一的约束线性化这是初学者最容易出错的地方。最初的直觉可能是如果 x_{ij}1则 d_{ij} R。但这不是一个线性约束。正确的线性化方法是引入一个足够大的常数M写成d_{ij} - R M * (1 - x_{ij})。当x_{ij}1时约束退化为d_{ij} - R 0即d_{ij} R当x_{ij}0时只要M足够大约束自动满足。而我之前使用的d_{ij} * x_{ij} R实际上要求的是“距离与0-1变量的乘积”小于等于R。当x_{ij}0时0 R恒成立当x_{ij}1时d_{ij} R。这同样是正确的线性化且更简洁。但需要注意的是这要求d_{ij}非负本题中距离显然满足。问题二的模型选择警车调度问题本质上是一个指派问题。我构建的模型是标准的二分图匹配模型。但需要注意的是原题是从20个平台中选13个去13个路口这等同于一个“不平衡指派问题”。我的约束设置每个路口必有一车、每个平台至多出一车、总车数等于13是正确的。也可以转化为一个标准的“平衡指派问题”即虚拟增加7个“虚拟路口”其到所有平台的距离为0这样就能用经典的匈牙利算法求解。用PuLP等通用求解器解0-1规划对于这个规模的问题13*20260个变量是瞬间完成的。6.3 求解器选择与性能调优我使用了PuLP默认的CBC求解器。对于这种规模的线性整数规划问题CBC完全够用。但如果问题规模扩大比如节点数上千可能会遇到求解速度慢的问题。这时可以尝试商业求解器PuLP支持Gurobi、CPLEX等。如果有许可证它们速度更快。调整求解参数例如设置时间限制、容忍间隙。模型简化对于问题一我们其实可以不建立庞大的0-1规划模型而是用聚类思想。因为目标是最小化最大距离这很像K-means聚类但距离是网络距离而非欧式距离。我们可以用启发式算法如反复将距离平台最远的节点分配给其他平台并调整可能快速得到一个近似最优解用于验证精确解的正确性。6.4 可视化中的信息过载与清晰表达最初的可视化尝试可能想把所有信息都塞进一张图网络、平台、管辖范围、调度路线。结果往往是一团乱麻。我的经验是分层展示先画基础网络和平台位置。然后用交互式工具如Plotly让用户可以点击选择某个平台高亮显示其管辖范围。最后再单独绘制调度路线图。简化边在绘制调度路线时直接连线如我示例中的虚线可能会误导因为它可能不是实际道路。更准确的做法是利用Floyd算法中记录的前驱节点信息回溯出最短路径的具体节点序列然后将这条路径上的所有边高亮绘制。这需要修改Floyd算法来记录路径信息计算量稍大但结果更精确。使用图布局算法原题给了坐标所以直接用。如果没有坐标需要用nx.spring_layout等算法自动计算节点位置但这样得到的布局可能与实际地理无关仅反映拓扑结构。复现经典赛题就像与过去的优秀选手隔空对话。用今天的工具重新求解不仅能确保答案的正确性更能深刻理解问题背后的数学模型和算法思想。从数据清洗、图论计算、数学建模到求解可视化这一整套流程正是解决当今许多现实世界优化问题如物流配送、网络规划、资源调度的缩影。把这个项目做扎实其意义远超一道题本身。
返回列表