Python 3.12 + Pandas 复现 MathorCup B题：共享单车时空分布与OD矩阵计算
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/Seaborn）：
import matplotlib.pyplot as plt import seaborn as sns plt.figure(figsize=(12, 6)) sns.lineplot(data=hourly_stats, x='hour', y='trip_count', marker='o') 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:00
2.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_value=0) # 标准化处理（出行比例） od_matrix_normalized = od_matrix.div(od_matrix.sum(axis=1), axis=0)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_value=0) 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_minutes=30): time_range = pd.date_range( start=clean_data['start_time'].min().floor('H'), end=clean_data['start_time'].max().ceil('H'), freq=f'{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[i+1]) ] od = window_data.groupby(['start_region', 'end_region']).size().unstack(fill_value=0) 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_size=0.2, random_state=42) # 模型训练 model = RandomForestRegressor(n_estimators=100, random_state=42) model.fit(X_train, y_train)4.2 调度路径优化
使用OR-Tools实现VRP：
from 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, labels=dict(x="Destination", y="Origin", color="Trip Ratio"), x=[f'Zone {i}' for i in range(1, 11)], y=[f'Zone {i}' for i in range(1, 11)], color_continuous_scale='Viridis') fig.update_layout(title='Normalized 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(lower=0) 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）
- 调度系统的容错机制（自动重试、故障转移）