数学建模与Python:共享单车系统建模-从简单模拟到参数优化

# 初始化状态：字典是最简单的数据结构state = {'campus': 100, 'metro': 20}

defbike_to_metro(state):"""将一辆车从校园移动到地铁站"""    state['campus'] -= 1    state['metro'] += 1

import numpy as npdefstep(state, p1, p2):"""模拟一个时间步（例如一分钟）"""# 模拟从校园出发if np.random.random() < p1:        bike_to_metro(state)# 模拟从地铁站出发if np.random.random() < p2:        bike_to_campus(state)

import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.animation import FuncAnimationimport matplotlib.patches as patches# ==========================================# 共享单车系统仿真模型# ==========================================# 1. 系统参数TOTAL_BIKES = 12P_CAMPUS_TO_METRO = 0.5# 校园 -> 地铁站 的概率 (每分钟)P_METRO_TO_CAMPUS = 0.4# 地铁站 -> 校园 的概率 (每分钟)DURATION = 60# 模拟时长 (分钟)# 2. 初始化状态# bikes_campus: 校园里的单车数# bikes_metro: 地铁站里的单车数state = {'campus': 10, 'metro': 2,'campus_history': [],'metro_history': [],'movements': [] # 记录每一帧的移动事件: None, 'c2m', 'm2c'}# 3. 模拟函数defrun_step(current_state):# 记录历史    current_state['campus_history'].append(current_state['campus'])    current_state['metro_history'].append(current_state['metro'])# 模拟移动    move = None# 检查从校园到地铁站的借车请求if np.random.random() < P_CAMPUS_TO_METRO:if current_state['campus'] > 0:            current_state['campus'] -= 1            current_state['metro'] += 1            move = 'c2m'# 检查从地铁站到校园的借车请求if np.random.random() < P_METRO_TO_CAMPUS:if current_state['metro'] > 0:            current_state['metro'] -= 1            current_state['campus'] += 1            move = 'm2c'    current_state['movements'].append(move)# 运行完整模拟以获取数据np.random.seed(42) # 固定种子保证结果可复现sim_data = {'campus': [10],'metro': [2],'movements': [None]}temp_state = {'campus': 10, 'metro': 2, 'campus_history': [], 'metro_history': [], 'movements': []}for _ in range(DURATION):    run_step(temp_state)sim_data['campus'].extend(temp_state['campus_history'])sim_data['metro'].extend(temp_state['metro_history'])sim_data['movements'].extend(temp_state['movements'])# ==========================================# 动态可视化 (Dashboard 风格)# ==========================================plt.style.use('dark_background')fig = plt.figure(figsize=(12, 7))gs = fig.add_gridspec(2, 2, height_ratios=[1.5, 1])# --- 上半部分：实时场景模拟 ---ax_map = fig.add_subplot(gs[0, :])ax_map.set_xlim(0, 100)ax_map.set_ylim(0, 40)ax_map.set_title("Real-time Bike Sharing Monitor", color='white', fontsize=16, pad=10)ax_map.axis('off')# 绘制地点# 校园 (左侧)campus_rect = patches.FancyBboxPatch((5, 5), 25, 30, boxstyle="round,pad=0.5",                                    fc='#27ae60', ec='#2ecc71', alpha=0.3)ax_map.add_patch(campus_rect)ax_map.text(17.5, 30, "CAMPUS", ha='center', color='#2ecc71', fontsize=12, fontweight='bold')text_campus_count = ax_map.text(17.5, 15, "", ha='center', color='white', fontsize=24, fontweight='bold')# 地铁站 (右侧)metro_rect = patches.FancyBboxPatch((70, 5), 25, 30, boxstyle="round,pad=0.5",                                   fc='#2980b9', ec='#3498db', alpha=0.3)ax_map.add_patch(metro_rect)ax_map.text(82.5, 30, "METRO STN", ha='center', color='#3498db', fontsize=12, fontweight='bold')text_metro_count = ax_map.text(82.5, 15, "", ha='center', color='white', fontsize=24, fontweight='bold')# 道路ax_map.plot([30, 70], [20, 20], '-', color='#7f8c8d', lw=4, zorder=0)ax_map.plot([30, 70], [20, 20], '--', color='white', lw=1, zorder=0)# 骑行者图标 (初始隐藏)rider = ax_map.text(50, 22, "🚴", ha='center', fontsize=30, color='#f1c40f', alpha=0)# --- 下半部分：数据图表 ---# 校园单车数量曲线ax_chart = fig.add_subplot(gs[1, :])ax_chart.set_xlim(0, DURATION)ax_chart.set_ylim(0, TOTAL_BIKES + 2)ax_chart.set_xlabel("Time (min)")ax_chart.set_ylabel("Available Bikes")ax_chart.grid(True, linestyle='--', alpha=0.2)line_campus, = ax_chart.plot([], [], '-', color='#2ecc71', lw=2, label='Campus')line_metro, = ax_chart.plot([], [], '-', color='#3498db', lw=2, label='Metro Station')ax_chart.legend(loc='upper right', frameon=False)# 时间指示器time_text = ax_chart.text(0.02, 0.9, "", transform=ax_chart.transAxes, color='white')defupdate(frame):# 1. 更新计数文本    c_count = sim_data['campus'][frame]    m_count = sim_data['metro'][frame]    text_campus_count.set_text(f"{c_count}")    text_metro_count.set_text(f"{m_count}")# 2. 更新骑行动画# 获取当前分钟发生的移动事件    move = sim_data['movements'][frame]if move == 'c2m': # 校园 -> 地铁        rider.set_text("🚴 ➤")        rider.set_x(40) # 简化的位置示意        rider.set_alpha(1)        rider.set_color('#2ecc71') # 绿色代表来自校园elif move == 'm2c': # 地铁 -> 校园        rider.set_text("◀ 🚴")        rider.set_x(60)        rider.set_alpha(1)        rider.set_color('#3498db') # 蓝色代表来自地铁else:        rider.set_alpha(0) # 没人骑车# 3. 更新曲线图    line_campus.set_data(range(frame+1), sim_data['campus'][:frame+1])    line_metro.set_data(range(frame+1), sim_data['metro'][:frame+1])    time_text.set_text(f"Time: {frame} min")return text_campus_count, text_metro_count, rider, line_campus, line_metro, time_text# 生成动画ani = FuncAnimation(fig, update, frames=DURATION+1, interval=100, blit=True)ani.save('../images/02_bike_share_dynamic.gif', writer='pillow', fps=10)print("Animation saved to ../images/02_bike_share_dynamic.gif")

动态演示效果：

• 上半部分：展示了校园（左）和地铁站（右）的实时场景。中间的骑行者图标表示正在发生的移动事件。
• 下半部分：展示了两个地点单车数量随时间的实时变化曲线。

4. 完善模型：处理“无车可用”的情况

上面的模型有一个明显的漏洞：如果校园已经没车了（Campus = 0），函数 bike_to_MetroStn 仍然会执行减一操作，导致车数变成负数。这在物理上是不可能的。

我们需要修改移动车辆的函数，加入检查逻辑。如果没车了，就记录一次“不满意的客户”（Unsatisfied Customer），直接返回，不再减车。

defbike_to_MetroStn(state):"""尝试将一辆车从校园移动到地铁站"""if state.Campus == 0:        state.Campus_empty += 1# 记录不满意的客户return    state.Campus -= 1    state.MetroStn += 1

这样，Campus_empty 和 MetroStn_empty 就成了我们评估系统服务质量的重要指标。

可视化展示与参数扫视

模拟结果可视化

运行 run_simulation(0.3, 0.2, 60)，我们可以得到一张图表，显示校园校园内单车数量随时间的变化。由于引入了随机性，每次运行生成的曲线都会有所不同，呈锯齿状波动。

参数扫视（Parameter Sweep）

仅仅运行一次模拟是不够的。我们想知道，不同的借车概率（p1）对系统服务质量到底有多大影响？比如，当校园的借车需求增加时，不满意的客户数量会如何变化？

这就需要进行参数扫视。我们让 p1 从 0 到 0.8 逐步变化，观察“校园不满意的客户数”这个指标。为了减少随机波动，我们在每个参数点上重复运行 20 次模拟取平均值。

你可以运行代码来生成参数扫视的动态演示，完整代码如下：

import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.animation import FuncAnimation# ==========================================# 参数扫视模拟# ==========================================TOTAL_BIKES = 12DURATION = 60NUM_SIMULATIONS = 20# 每个参数值重复模拟的次数 (用于平滑曲线)P2 = 0.4# 地铁站 -> 校园 的固定概率# 借车概率 (校园 -> 地铁) 范围：0 到 0.8p1_array = np.linspace(0, 0.8, 41)# 存储结果：每个 p1 对应的不满意客户数 (平均值)sweep_results = []current_p1_index = 0defrun_simulation(p1, p2, num_steps):    state = {'campus': 10, 'metro': 2, 'unhappy': 0}for _ in range(num_steps):# 校园 -> 地铁if np.random.random() < p1:if state['campus'] > 0:                state['campus'] -= 1                state['metro'] += 1else:                state['unhappy'] += 1# 地铁 -> 校园if np.random.random() < p2:if state['metro'] > 0:                state['metro'] -= 1                state['campus'] += 1# 这里简化，不统计地铁站的不满意客户return state['unhappy']# ==========================================# 动态可视化# ==========================================plt.style.use('dark_background')fig, ax = plt.subplots(figsize=(10, 6))ax.set_xlim(0, 0.8)ax.set_ylim(0, 15)ax.set_xlabel("Request Rate at Campus (p1)")ax.set_ylabel("Average Unhappy Customers")ax.set_title("Impact of Request Rate on Service Quality", fontsize=14, color='white', pad=20)ax.grid(True, linestyle='--', alpha=0.3)# 绘制曲线line, = ax.plot([], [], 'o-', color='#e74c3c', lw=2, markersize=6, label='Unhappy Customers')ax.legend(loc='upper left', frameon=False)# 动态文本text_param = ax.text(0.05, 0.8, "", transform=ax.transAxes, color='white', fontsize=12)text_result = ax.text(0.05, 0.75, "", transform=ax.transAxes, color='#e74c3c', fontsize=14, fontweight='bold')# 进度条背景ax_progress = fig.add_axes([0.15, 0.02, 0.7, 0.03]) # [left, bottom, width, height]ax_progress.set_xlim(0, 1)ax_progress.set_ylim(0, 1)ax_progress.axis('off')rect_bg = plt.Rectangle((0, 0), 1, 1, fc='#34495e', ec=None)ax_progress.add_patch(rect_bg)rect_progress = plt.Rectangle((0, 0), 0, 1, fc='#2ecc71', ec=None)ax_progress.add_patch(rect_progress)ax_progress.text(0.5, 0.5, "Simulation Progress", ha='center', va='center', color='white', fontsize=8)defupdate(frame):global current_p1_indexif current_p1_index < len(p1_array):        p1 = p1_array[current_p1_index]# 运行多次模拟取平均值        total_unhappy = 0for _ in range(NUM_SIMULATIONS):            total_unhappy += run_simulation(p1, P2, DURATION)        avg_unhappy = total_unhappy / NUM_SIMULATIONS        sweep_results.append(avg_unhappy)# 更新曲线        line.set_data(p1_array[:len(sweep_results)], sweep_results)# 更新文本        text_param.set_text(f"Current p1 (Campus->Metro): {p1:.2f}")        text_result.set_text(f"Unhappy Customers: {avg_unhappy:.1f}")# 更新进度条        progress = (current_p1_index + 1) / len(p1_array)        rect_progress.set_width(progress)        current_p1_index += 1return line, text_param, text_result, rect_progress# 生成动画ani = FuncAnimation(fig, update, frames=len(p1_array)+10, interval=100, blit=True)ani.save('../images/02_bike_share_sweep.gif', writer='pillow', fps=10)print("Animation saved to ../images/02_bike_share_sweep.gif")

动态参数扫视效果：

通过这张动态图，我们可以清晰地看到：

1. 当 p1 较小时（需求低），不满意客户数维持在 0 附近。
2. 随着 p1 增大，曲线开始出现非线性上升。这意味着一旦需求超过某个临界点，服务质量会急剧下降。
3. 这种“相变”现象提示我们，在设计系统容量时，必须留有足够的余量来应对高峰需求。

总结

通过这个共享单车系统的案例，我们完整地体验了建模的过程：

1. 从简单开始：先建立一个能跑通的最小模型（State 和基本移动函数）。
2. 迭代改进：发现问题（负数车辆），修改代码（增加条件判断），引入度量指标（不满意客户）。
3. 探索分析：使用参数扫视来分析系统行为随参数变化的规律。

这种迭代式建模（Iterative Modeling） 的方法，不仅适用于简单的共享单车系统，也是解决复杂工程和科学问题的有效路径。