Python 3.12 Pandas 复现 MathorCup B题共享单车时空分布与OD矩阵计算实战指南1. 项目背景与数据准备共享单车作为城市短途出行的重要工具其运营效率直接影响用户体验和企业效益。2017年MathorCup数学建模B题以共享单车调度优化为背景要求参赛者分析时空分布特征并构建OD矩阵。本文将使用Python 3.12和Pandas生态系统完整复现该案例。核心数据集应包含以下字段start_time: 骑行开始时间datetime格式end_time: 骑行结束时间datetime格式start_region: 出发区域编码整数end_region: 到达区域编码整数user_id: 用户唯一标识符典型数据加载与预处理代码import pandas as pd import numpy as np # 加载原始数据 raw_data pd.read_csv(bike_trips.csv, parse_dates[start_time, end_time], dtype{start_region: int8, end_region: int8}) # 计算骑行时长分钟 raw_data[duration] (raw_data[end_time] - raw_data[start_time]).dt.total_seconds() / 60 # 过滤异常数据 clean_data raw_data[ (raw_data[duration] 1) (raw_data[duration] 180) (raw_data[start_region].between(1, 10)) (raw_data[end_region].between(1, 10)) ].copy()提示实际项目中建议使用pandas-profiling或ydata-profiling生成数据质量报告快速识别异常值和缺失值。2. 时空分布特征分析2.1 时间维度分析小时级骑行量统计# 提取小时特征 clean_data[hour] clean_data[start_time].dt.hour # 分时段统计 hourly_stats clean_data.groupby(hour).agg( trip_count(user_id, count), avg_duration(duration, mean) ).reset_index()可视化时段分布需配合Matplotlib/Seabornimport matplotlib.pyplot as plt import seaborn as sns plt.figure(figsize(12, 6)) sns.lineplot(datahourly_stats, xhour, ytrip_count, markero) plt.title(Hourly Bike Trip Distribution) plt.xlabel(Hour of Day) plt.ylabel(Trip Count) plt.grid(True) plt.show()典型时间分布特征早高峰7:00-9:00通勤需求午间小高峰11:00-13:00午餐出行晚高峰17:00-19:00下班通勤夜间低谷22:00-次日5:002.2 空间维度分析区域热度排名departure_rank clean_data[start_region].value_counts().sort_index() arrival_rank clean_data[end_region].value_counts().sort_index() region_stats pd.DataFrame({ departures: departure_rank, arrivals: arrival_rank, net_flow: arrival_rank - departure_rank })区域类型特征描述净流入区商业中心、交通枢纽净流出区住宅区、学校平衡区混合功能区3. OD矩阵构建与优化3.1 基础OD矩阵计算# 构建原始OD矩阵 od_matrix clean_data.groupby([start_region, end_region]).size().unstack(fill_value0) # 标准化处理出行比例 od_matrix_normalized od_matrix.div(od_matrix.sum(axis1), axis0)3.2 时空增强OD矩阵分时段OD矩阵示例def build_time_slice_od(df, time_ranges): results [] for name, (start_h, end_h) in time_ranges.items(): slice_data df[df[hour].between(start_h, end_h)] od slice_data.groupby([start_region, end_region]).size().unstack(fill_value0) od[time_slice] name results.append(od) return pd.concat(results) time_ranges { morning: (7, 9), noon: (11, 13), evening: (17, 19) } enhanced_od build_time_slice_od(clean_data, time_ranges)3.3 动态OD矩阵优化基于滑动窗口的实时更新def dynamic_od_update(window_minutes30): time_range pd.date_range( startclean_data[start_time].min().floor(H), endclean_data[start_time].max().ceil(H), freqf{window_minutes}T ) dynamic_ods [] for i in range(len(time_range)-1): window_data clean_data[ (clean_data[start_time] time_range[i]) (clean_data[start_time] time_range[i1]) ] od window_data.groupby([start_region, end_region]).size().unstack(fill_value0) od[window_start] time_range[i] dynamic_ods.append(od) return pd.concat(dynamic_ods)4. 调度优化模型实现4.1 需求预测模型from sklearn.ensemble import RandomForestRegressor from sklearn.model_selection import train_test_split # 特征工程 features pd.get_dummies(clean_data[[hour, start_region, end_region]]) target clean_data[duration] # 训练测试分割 X_train, X_test, y_train, y_test train_test_split( features, target, test_size0.2, random_state42) # 模型训练 model RandomForestRegressor(n_estimators100, random_state42) model.fit(X_train, y_train)4.2 调度路径优化使用OR-Tools实现VRPfrom ortools.constraint_solver import routing_enums_pb2 from ortools.constraint_solver import pywrapcp def optimize_route(demand_matrix, distance_matrix): # 创建距离回调函数 def distance_callback(from_index, to_index): from_node manager.IndexToNode(from_index) to_node manager.IndexToNode(to_index) return distance_matrix[from_node][to_node] # 初始化路由模型 manager pywrapcp.RoutingIndexManager( len(distance_matrix), 1, 0) # 1辆车0为仓库 routing pywrapcp.RoutingModel(manager) transit_callback_index routing.RegisterTransitCallback(distance_callback) routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index) # 设置搜索参数 search_parameters pywrapcp.DefaultRoutingSearchParameters() search_parameters.first_solution_strategy ( routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC) # 求解 solution routing.SolveWithParameters(search_parameters) # 提取最优路径 index routing.Start(0) route [] while not routing.IsEnd(index): route.append(manager.IndexToNode(index)) index solution.Value(routing.NextVar(index)) route.append(manager.IndexToNode(index)) return route5. 可视化与结果分析5.1 热力图可视化import plotly.express as px fig px.imshow(od_matrix_normalized, labelsdict(xDestination, yOrigin, colorTrip Ratio), x[fZone {i} for i in range(1, 11)], y[fZone {i} for i in range(1, 11)], color_continuous_scaleViridis) fig.update_layout(titleNormalized OD Matrix) fig.show()5.2 动态调度模拟def simulate_dynamic_scheduling(od_series, initial_inventory): inventory initial_inventory.copy() records [] for ts, od in od_series.items(): # 记录当前库存 record inventory.copy() record[timestamp] ts records.append(record) # 模拟骑行消耗 for origin in od.index: for dest in od.columns: trips od.loc[origin, dest] inventory[origin] - trips inventory[dest] trips # 库存修正防止负值 inventory inventory.clip(lower0) return pd.DataFrame(records).set_index(timestamp)6. 工程化实践建议性能优化技巧对大规模数据使用pandas.DataFrame.eval()进行链式运算使用numba加速数值计算密集型函数对分类数据使用category数据类型减少内存占用部署架构graph TD A[实时数据流] -- B(Kafka消息队列) B -- C{Spark流处理} C -- D[OD矩阵计算] C -- E[需求预测] D -- F[调度优化引擎] E -- F F -- G[调度指令下发]异常处理模式class BikeDataException(Exception): pass def validate_od_matrix(matrix): if matrix.isnull().any().any(): raise BikeDataException(OD矩阵包含空值) if (matrix 0).any().any(): raise BikeDataException(OD矩阵包含负值) return True在实际项目中我们还需要考虑与GIS系统的集成使用GeoPandas处理空间数据实时数据管道建设Apache Kafka Spark Streaming调度系统的容错机制自动重试、故障转移
